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AN  ADAPTIVE  FINITE  DIFFERENCE  METHOD  FOR 
HYPERBOLIC  SYSTEMS  IN  ONE  SPACE  DIMENSION 


Abstract  Many  problems  of  physical  interest  have  solutions  which  are  gen¬ 
erally  quite  smooth  in  a  large  portion  of  the  region  of  interest,  but  have  local 
phenomena  such  as  shocks,  discontinuities  or  large  gradients  which  require 
much  more  accurate  approximations  or  finer  grids  for  reasonable  accuracy. 
Examples  are  atmospheric  fronts,  ocean  currents,  and  geological  discontinui¬ 
ties. 

>>In  this  paper  we  develop  and  partially  analyze  an  adaptive  finite  difference 
mesh  refinement  algorithm  for  the  initial  boundary  value  problem  for  hypdfbalic 
systems  in  one  space  dimension.  The  method  uses  clusters  of  uniform  grids 
which  can  "move"  along  with  pulses  or  steep  gradients  appearing  in  the  calcula¬ 
tion.  and  which  are  superimposed  over  a  uniform  coarse  grid.  Such  refinements 
are  created,  destroyed,  merged,  separated,  recursively  nested  or  moved  based 
on  estimates  of  the  local  truncation  error.  We  use  a  four-way  linkeu  tree  and 
sequentially  allocated  deques  (double-ended  queues)  to  perform  these  opera¬ 
tions  efficiently.  The  local  truncation  error  is  estimated  using  a  three-step 
Richardson  extrapolation  procedure  in  the  interior  of  the  region,  and  differences 
at  the  boundaries.  Our  algorithm  was  implemented  using  a  portable,  extensible 
Fortran  preprocessor,  to  which  we  added  records  and  pointers. 

The  method  is  applied  to  two  model  problems:  the  second  order  wave  equa¬ 
tion  with  counterstreaming  Gaussian  pulses,  and  the  Riemann  shock-tube  prob¬ 
lem.  For  both  problems  our  algorithm  is  shown  to  be  three  to  five  times  more 
efficient  (in  computing  time)  than  the  use  of  a  uniform  coarse  mesh,  for  the 
same  accuracy.  Furthermore,  to  our  knowledge,  our  algorithm  is  the  only  one 
which 
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1.  Introduction.  Many  problems  of  physical  interest  have  solutions  which 

are  smooth  in  a  large  portion  of  the  region  of  interest,  but  have  local 

phenomena  such  as  shocks,  discontinuities  or  large  gradients  which  require 

much  more  accurate  approximation  or  finer  meshes  for  reasonable  accuracy. 

* 

Examples  of  this  are  atmospheric  fronts,  ocean  currents,  geological  discontinui¬ 
ties.  and  storm  surges. 

When  the  positions  of  the  gradients  are  known  a  priori,  and  are  independent 
of  time,  one  can  use  coordinate  transformations,  a  technique  used  extensively  in 
aerodynamic  computations,  e.g.,  Steger  and  Chaussee  [44].  If  the  position  of 
the  gradients  changes  as  an  a  priori  function  of  time,  the  mapping  function  can 
change  with  time  (e.g. ,  flow  past  a  helicopter  blade).  However,  when  the  manner 
in  which  the  gradients  move  is  not  known  in  advance,  this  technique  cannot  be 
used. 

The  use  of  a  fine  mesh  throughout  the  entire  calculation  region  usually 
requires  too  much  computer  time  and/or  storage.  An  alternative  method  is  to 
use  an  underlying  coarse  mesh  for  the  entire  region,  and  to  superimpose  a  fine 
grid,  or  grids,  on  the  region(s)  where  the  solution  is  varying  rapidly.  The  crucial 
difficulty  is  that  the  refined  region(s)  must  then  move  along  with  the  rapidly 
varying  portion  of  the  solution,  at  all  times  enclosing  this  portion. 

The  necessity  for  this  was  illustrated  by  Browning,  Kreiss  and  Oliger  [7]. 
They  computed  the  numerical  solution  of  an  initial  boundary  value  problem 
whose  exact  solution  was  a  rapidly  oscillating  sine  wave.  The  wave  was  accu¬ 
rately  represented  on  the  fine  part  of  the  spatial  grid,  but  when  it  passed  into 
the  coarse  mesh,  it  was  almost  completely  obliterated.  From  this  it  is  clear  that 
the  rapidly  varying  part  of  the  solution  must  not  be  allowed  to  escape  the 
refinement  region. 
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Adaptive  methods  have  been  used  in  other  areas  of  numerical  analysis  for 
many  years;  some  references  are  given  in  [5].  The  study  of  adaptive  algorithms 
for  time-dependent  partial  differential  equations  is  a  very  active  research  area. 
We  will  now  list  several  properties  which  distinguish  our  method  from  others, 
and  at  the  same  time  list  a  few  of  the  many  references  to  other  work. 

Wc  should  first  mention  that  the  methods  of  Gropp  [18],  Berger  [2]  and  the 
author  all  follow  the  strategy  of  Budnik  and  Oliger  [8],  Oliger  [33],  and  Berger, 
Gropp  and  Oliger  [3],  Gropp's  paper  is  a  test  of  the  feasibility  of  this  method  in 
two  dimensions.  Berger's  method  is  for  hyperbolic  problems  in  one  and  two 
space  dimensions.  It  uses  the  same  method  as  ours  for  placing  refinements  (a 
method  which  is  more  general  than  Gropp's)  and,  like  ours,  allows  fine  grids  to 
be  recursively  nested.  Berger's  method  allows  refinements  to  have  arbitrary 
orientation  in  two  dimensions.  Her  method  does  not  treat  boundary  conditions 
as  ours  does.  Our  data  structure  is  somewhat  different  than  hers,  but  cannot  be 
generalized  to  more  than  one  space  dimension. 

The  most  important  difference  between  our  algorithm  and  most  others  is  in 
efficiency.  On  model  problems  our  algorithm  is  three  to  five  times  more 
efficient  than  using  a  uniform  mesh  which  produces  the  same  accuracy.  (Berger 
obtains  similar  efficiency  factors,  while  Gropp’s  were  somewhat  less.)  Most  other 
methods  fail  to  discuss  this  matter,  and  fail  to  give  computer  times  for  their 
methods.  An  exception  is  Hu  and  Schiesser  [24],  who  found  their  method  to  be 
approximately  1.5  times  as  efficient  on  a  test  problem. 

The  second  distinguishing  feature  of  our  method  is  that  it  is  designed  for 
hyperbolic  systems.  Many  other  methods  claim  to  work  on  both  parabolic  and 
hyperbolic  problems,  but  their  examples  are  only  given  for  the  former.  In  con¬ 
trast,  our  method  can  be  used  for  parabolic  problems,  but  it  is  not  yet  efficient 
for  them  since  it  uses  explicit  time  steps.  Our  method  does  not  require  us  to 
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know  the  position  or  the  direction  of  characteristics. 

We  can  classify  adaptive  algorithms  as  either  (a)  moving  (or  deformed)  grid 
methods,  or  (b)  local  refinement  methods.  In  the  former,  a  fixed  number  of  grid 
points  are  present,  and  they  move  with  the  solution  (possibly  via  coordinate 
transformations).  In  more  than  one  space  dimension  this  can  lead  to  distorted 
grids.  The  other  approach  uses  a  stationary  coarse  grid,  on  which  are  overlayed 
or  interpolated  finer  grid  points.  To  follow  a  wave,  grid  points  are  created  or 
destroyed  as  necessary,  but  grid  points  are  not  actually  moved.  Our  method  is  a 
local  refinement  method,  and  our  refined  grids  can  be  created,  destroyed, 
merged,  separated  or  recursively  nested.  In  some  problems  there  is  a  shock  or 
gradient  that  needs  to  be  resolved  at  one  time,  and  several  at  another  time  (for 
example,  see  our  Riemann  shock  tube  problem  in  Section  B).  A  method  with  a 
fixed  number  of  points  may  use  too  many  or  too  few  part  of  the  time.  If  too  few 
are  used,  the  grid  points  may  also  undergo  sudden  transitions  when  another 
steep  gradient  appears.  Hedstrom  and  Rodrigue  [23]  discuss  other  advantages 
and  disadvantages  of  these  approaches.  The  methods  of  Davis  and  Flaherty  [11], 
Miller  [17],  [32]  Dwyer.  Kee  and  Sanders  [13],  Rai  and  Anderson  [39],  Brackbill 
and  Saltzman  [8],  and  White  [47]  are  moving  grid  methods.  The  methods  of 
Berger,  Gropp,  Gannon  [16],  Hu  and  Schiesser  Lam  and  Simpson  [30],  and  the 
author  are  local  refinement  methods. 

Among  the  local  refinement  methods  we  can  distinguish  those  in  which  fine 
grid  points  are  interpolated,  or  patched,  into  the  coarse  grid,  from  those  in 
which  the  fine  grid  points  overlay  (and  are  separate  from)  the  coarse  grid.  The 
methods  of  Lam  and  Simpson  and  Hu  and  Schiesser  are  examples  of  the  former, 
while  the  methods  of  Gropp,  Berger  and  the  author  all  use  the  latter  approach. 
An  advantage  of  the  overlay  approach  is  its  adaptability  to  parallel  or  pipeline 
processors.  The  solution  can  be  advanced  (and  the  error  estimated)  almost 
independently  on  different  grids. 


The  fourth  distinguishing  feature  of  our  method  is  that  it  refines  the  mesh 
in  both  space  and  time.  Most  other  adaptive  methods  use  a  time  step  which  is 
the  same  at  every  spatial  mesh  point  at  a  given  time.  Since  we  use  explicit 
methods,  it  is  essential  to  allow  finer  time  steps  in  refinements  than  in  the 
coarse  mesh.  Otherwise,  stability  would  force  us  to  use  tiny  time  steps  over  the 
whole  spatial  region.  To  our  knowledge,  only  the  methods  of  the  author,  Berger, 
and  Gropp  allow  this. 

Closely  related  to  this  is  the  fifth  distinguishing  feature  of  our  method:  the 
use  of  nontraditional  (for  numerical  methods)  data  structures.  This  is  neces¬ 
sary  to  efficiently  implement  variable  numbers  of  refined  grids.  Berger  has 
implemented  more  complicated  data  structures  in  two  space  dimensions.  Other 
adaptive  methods  using  nontrivial  data  structures  were  developed  by  Gannon 
and  by  Rheinboldt  and  Mesztenyi  [40]. 

A  sixth  notable  feature  of  our  method  is  our  criterion  for  mesh  placement. 
We  arrange  our  mesh  points  to  "approximately  equidistribute"  the  local  trunca¬ 
tion  error  (de  Boor  [12]  and  Pereyra  and  Sewell  [37]).  Other  methods  use  meas¬ 
ures  such  as  the  gradient,  arclength,  or  second  derivative.  We  feel  that  these  ad 
hoc  methods  may  work  in  certain  instances,  but  are  not  sufficiently  general. 
Other  methods  that  use  the  local  truncation  error  are  those  of  Pierson  and 
Kutler  [30],  Rai  and  Anderson,  and  Berger.  Gannon  uses  an  analogous  criterion 
[l]  for  finite  elements. 

A  unique  feature  of  our  method  is  the  adaptive  treatment  of  time- 
dependent  boundary  conditions.  Implementing  this  requires  local  error  esti¬ 
mates  at  the  boundary,  a  variable  number  of  grid  points,  and  nontrivial  data 
structures.  This  is  obviously  important  in  applications  like  limited  area  weather 
forecasting. 
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The  final  distinguishing  feature  of  our  algorithm  is  its  method  of  implemen¬ 
tation.  Fortran  is  inconvenient  because  it  lacks  data  structures  and  control 
structures,  and  other  languages  have  portability  problems.  So  we  used  an 
extensible,  portable  Fortran  preprocessor,  to  which  we  added  pointers  and 
records.  Stein  [45]  and  Gropp  [  19]  have  also  considered  languages  for  adaptive 
algorithms. 

We  now  summarize  what  is  contained  in  the  rest  of  this  paper.  In  Section  2 
we  describe  our  adaptive  mesh-refinement  algorithm  in  detail.  We  first  describe 
the  continuous  problem,  the  usual  first  order  hyperbolic  system  on  a  strip  in 
one  space  dimension.  Next  we  describe  our  mesh  structure.  A  detailed  descrip¬ 
tion  of  the  algorithm  is  then  provided,  including  techniques  at  boundaries  and  at 
interfaces  between  coarse  and  fine  meshes.  Section  3  discusses  the  estimation 
of  the  local  truncation  error  in  the  interior  of  refinements,  at  coarse /fine  inter¬ 
faces.  and  at  boundaries. 

In  Section  4  we  give  a  brief  discussion  of  stability.  Section  5  summarizes 
partial  results  in  [5],  which  show  how  mesh  refinement  affects  the  rate  of  con¬ 
vergence.  Section  6  describes  the  data  structures  we  used  to  implement  the 
algorithm.  Section  7  discusses  the  Fortran  preprocessor  language  we  used.  Sec¬ 
tion  8  provides  computational  results.  As  model  problems  we  used  the  first 
order  wave  equation,  the  second  order  wave  equation  with  counterstreaming 
Gaussian  pulses,  and  the  Riemann  shock  tube  problem. 

2.  Uesh  Structure  and  Solution  Algorithm.  In  this  section  we  describe  the 
differential  equations  to  be  approximated.  Then  we  give  a  recursive  description 
of  our  mesh  structure  and  of  our  algorithm  for  advancing  the  solution. 

2.1.  The  Continuous  Problem.  Let  G  denote  the  spatial  interval  a  <  x  «  6 . 
For  purposes  of  analysis  we  will  consider  a  linear  first  order,  one  (space)- 


6 


dimensional,  n  x  n  hyperbolic  system 

Lu  -  Ui  -  A{z,t)ug  -  B(x,t)u  -  F{x,t),  (2.1) 

on  a  "vertical"  strip  Cl  x  \t  fe  0{.  with  initial  condition 

u(x,0 )  =  /(*),  x  e  Cl,  (2.2) 

and  boundary  conditions 

u^o.t)  =  5n(i)un(a,0  +  Pi(0.  *  ^  0,  (2.3) 

un(b,0  =  ^i(0«I(*>.0  +  ^2(0.  0.  (2.4) 


Here  A  and  B  are  n  x  n  matrices  and  F  is  an  n-vector.  We  have,  as  usual, 
assumed  that  A  has  already  been  transformed  into  diagonal  form  by  a  nonsingu¬ 
lar  uniformly  bounded  similarity  transformation  T(x,t),  so  that 


with  T{x.t)  and  T{x,t)~l  uniformly  bounded,  and 

A1  =  diagta,  «2 . */)  <  0. 

An  =  diag(*v+i,  K/+2 . Kn)  >  0, 

U 1  =  (u1(x,f),  UZ(x,t ) . Uj{x,t))T, 

Ua  =  (uy  +  1(x,f).  Uj+2(x,t) . Unix.t))7, 

and 

5,  e  Su  e 

By  far  the  most  important  restriction  is  that  our  problem  has  only  one 
space  dimension.  Even  in  two  space  dimensions  the  problem  has  severe  addi¬ 
tional  difficulties,  such  as  irregular  geometries,  orientation  of  refinements,  pat¬ 
tern  recognition,  the  need  for  more  complicated  data  structures,  and  topology 
of  boundaries.  The  restriction  to  hyperbolic  behavior  insures  that  we  can  use 
explicit  time  steps.  This  assumption  greatly  simplifies  both  the  error  estimation 
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and  the  manipulation  of  moving  meshes.  However,  many  computational  prob¬ 
lems  in  fluid  dynamics  and  elsewhere  are  of  this  type. 

We  have  assumed  that  the  matrix  A  is  in  diagonal  form  because  this  makes 
it  easier  to  write  down  boundary  conditions  (2.3)-(2.4)  which  yield  a  well-posed 
problem.  We  have  assumed  the  problem  is  linear  so  that  we  can  analyze  local 
error  estimation  and  convergence  (in  fact,  we  shall  assume  constant  coefficients 
in  Sections  3  and  5).  Neither  of  these  assumptions  is  necessary  in  practice,  as 
shown  by  computations  in  Section  B. 

We  assume  that  the  overall  phenomena  being  studied  are  such  that,  except 
for  relatively  small  regions,  a  coarse  uniform  mesh  is  sufficient  to  resolve  them. 
We  further  assume  these  small  regions  change  with  time  in  a  way  which  cannot 
conveniently  be  determined  a  priori.  We  also  assume  that,  at  any  time,  these 
small  regions  are  approximately  the  same  for  all  solution  components.  In  other 
words,  if  the  differential  equations  describe  velocity  and  pressure,  then  large 
pressure  gradients  occur  in  approximately  the  same  regions  where  large  velo¬ 
city  gradients  occur.  Thus  we  can  use  the  same  refinements  for  each  com¬ 
ponent  of  the  solution  vector.  (The  assumptions  in  this  paragraph  are  necessary 
only  for  efficiency.  The  method  will  work  without  them,  but  it  might  refine  too 
large  a  portion  of  the  region  ) 

For  purpose  of  analysis  we  assume  that  the  solution  is  smooth.  This  means 
that  there  are  no  corner  discontinuities,  i.e.,  it1  (a,0)  =  fl  (a)  and 
un(6,0)  =  /n(6).  Furthermore,  it  means  there  are  no  shocks  present.  These 
assumptions  enable  us  to  estimate  the  local  truncation  error  using  asymptotic 
expansions.  These  restrictions  are  not  always  necessary  in  practice  (see  prob¬ 
lem  P3  in  Section  8).  but  then  our  method  is  not  supported  by  analysis. 

2.2.  Mesh  Structure.  We  will  now  formally  describe  otir  refined  grids.  We 
will  compute  on  a  basic  rectangle  R  =  0  x  [  T 1  wh‘  is  subdivided  into 
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subrectangles  on  each  of  which  is  imposed  a  uniform  grid  (see  Fig.  1). 

We  will  proceed  recursively  by  "levels  of  refinement".  The  word  level  in  this 
context  refers  not  to  the  time  level,  but  to  how  fine  a  grid  spacing  is.  Finer  grids 
will  have  higher  levels.  Let  A  be  the  maximum  level.  On  each  level 

£=0,1 . A-l  we  will  introduce  a  finite  number  of  space-time  refinement 

rectangles  or  boxes  Blv  contained  in  rectangle  R.  (All  such  rectangles  will  be 
solid,  that  is,  they  include  both  interior  and  boundary.)  Each  such  rectangle  will 
have  sides  parallel  to  the  coordinate  axes,  and  for  £  i  1  each  £-th  level  rectangle 
must  lie  entirely  in  an  £-l-st  level  rectangle.  Furthermore,  no  two  £-th  level 
rectangles  can  overlap.  The  boundary  of  each  £-th  level  rectangle  will  be  the 
boundary  of  a  uniform  £+l-st  level  (space-time)  grid.  All  £  +  l-st  level  grids  will 
have  the  same  space  and  time  steps.  Loosely  speaking,  an  f  +  l-st  level 
refinement  is  one  of  these  grids  viewed  at  a  fixed  time. 

To  prime  the  recursive  pump,  wc  will  define  the  zero-th  level  spatial  divi¬ 
sion  points  of  the  interval  [a,  b]  as  the  sequence  of  points  <x§  =  a,  x?  =  b>. 
Similarly,  the  zero-th  level  time  division  points  of  the  interval  [0,  T]  comprise 
the  sequence  .1°  -  0.  £?  =  T>.  Let  h0  =  b  -  a  and  kQ=  T  be  the  zero-th  level 
space  and  time  steps,  respectively.  We  define  U°  as  the  set  of  four  corner  points 
of  the  rectangle  R. 

For  £  =  0,  1 . A-l  we  now  form  the  l-th  level  partition  Pi  of  [0,  7*],  which 

is  a  subsequence  of  the  time  division  points  <fm>' 

0  =  £o  <  t}nj  <  t}nz  <  •  •  <  £mfj_ J  <  =  T.  (2.5) 

We  have  omitted  from  the  notation  the  dependence  of  the  subsequence  <mi>  on 
£.  For  £  =  0,  F0  is  identical  to  the  sequence  <£m>  of  time  division  points.  Thus 
s0  =  1  and  m^O)  =  1.  For  £  1,  must  contain  as  a  subsequence  the  points  in 

the  partition  Pt.r 
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This  partition  divides  the  region  R  into  l-th  level  horizontal  strips 

Si  i  =  1,  2 . sj.  For  1=0  the  only  such  strip  5,°  is  identical  to  the  rectangle 

R.  For  l  >  1  each  of  these  strips  is  contained  in  an  i-l-st  level  strip,  since  Pi-i 
is  a  subsequence  of  Pt .  The  partition  points  are  the  times  when  we  adjust  the 
mesh.  They  must  be  chosen  before  the  solution  of  the  problem.  (The  partitions 
and  strips  for  l  >  1  could  be  dispensed  with  if  we  never  adjust  the  mesh  between 
coarse  time  steps.)  For  l  =  1  we  have  shown  in  Fig.  1  the  time  division  points 
<t^>  and  the  partition  Pi  with  m( }  =  0,  =  2,  =  3.  and  m.j  =  5.  We  shall 

denote  the  partition  points  Pi  alternatively  as 

0  =  t°  <<»  <  tz  <  ■  ■  •  <  f*  =  T,  (2.6) 

and  denote  the  first  level  strips  S by  St.  Thus  Si  =  fl  x  [ft_1,  f‘]  for 
i  =  1.  ■  ■  ,s. 

We  will  now  introduce  a  set  of  zero  or  more  nonoverlapping  l-th  level  (solid) 
refinement  rectangles 

\Blv,  v=  1.  2 . gil 

(If  <71  =  0  the  recursion  ends.)  There  is  only  one  zero-th  level  rectangle  Bf ,  and 
it  is  identical  to  the  rectangle  R.  For  l  at  1,  each  rectangle  B\,  is  required  to  lie 
entirely  in  some  l -1-st  level  refinement  rectangle  B\rl  •  The  latter  will  be  called 
a  parent  of  the  former.  Each  such  rectangle  Blv  will  have  horizontal  sides  whose 
f-coordinates  are  required  to  be  adjacent  members  of  the  partition  Pt  (2.5). 
That  is,  the  horizontal  sides  of  the  rectangle  are  the  same  as  the  the  horizontal 
sides  of  the  l-th  level  strip  in  which  it  is  contained.  Since  Pi- j  is  a  subsequence 
of  Pi  for  l  1,  we  are  guaranteed  that  B\,  is  "vertically  contained"  in  its  parent. 

For  l  1,  the  x-coordinates  of  the  vertical  sides  of  rectangle  Blv  can  be  any 
l-th  level  spatial  division  point,  so  long  as  Blv  is  "horizontally  contained"  in  its 
parent  In  other  words,  let  its  parent  B have  left  and  right  vertical  sides  with 
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coordinates 

z  =  Zjjii)  =  *^t_1)a(ir)  =  a  +  tn-Mn) 

and 

*  =  *l(»>  =  =  »  +  ^-ir¬ 

respectively.  (Here  a(n)  and  cj(tt)  are  nonnegative  integers.)  Then  for  the  coor¬ 
dinates  xla(v)  and  2[,(y)  of  the  left  and  right  vertical  sides  of  rectangle  £?{,,  we 
require 

xa(n)  ^  xa{v)  ^  xu[v)  ^  31  L(it)  ■  (2-7a) 

i.e., 

r)  <  a(u)  <  o(v)  <  N^l~^u(n).  (2.7b) 

If  the  left  (right)  edge  of  the  parent  rectangle  lies  on  the  left  (right)  edge  of  the 
region  R,  tnen  we  allow  the  leftmost  (rightmost)  inequality  to  become 

For  i&0,  let  7V(I)  and  (the  i-th  level  spatial  and  time 

refinement  ratios )  be  integers  greater  than  one.  (Fig.  1  shows  the  case  =  3 
and  =  2  for  l  >  1.)  Let  hl+l  =  \/ and  fe(+1  =  kt /  be  the  l+l-st  level 
space  and  time  steps,  respectively.  We  now  define  the  sequences  of  (uniform) 
£+  J-st  level  spatial  and  time  division  points 

<xj+1  =  a  +  ;7i{  + j  -  0,1 . 

M=0 

and 

T1*1  =  <C  =  mJb,+1:  m  =  0,1 . 

M=0 

of  the  intervals  Q  and  [0,  7’J,  respectively.  They  are  respectively  N ^  and 
times  as  fine  as  the  £-th  level  ones.  The  set  of  all  points 


ul"  =  K Xj",  **♦«)! 

occupies  the  entire  rectangle  R  =  Q  x  [0,  T].  The  subset  of  these  points  con- 
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tained  in  the  (solid)  refinement  rectangle  Blv  is  defined  to  be  the  (f  +  l)-st  Level 
(space-time)  grid  G\fl  occupying  Blv.  More  specifically,  if  B l  occupies  the  i-th 
level  horizontal  strip  Si,  then  Gf,*1  consists  of  that  subset  of  U1*1  whose  *  com¬ 
ponents  have  subscripts 

j  =  a(v)Nw,  a(v)N^+l . u(v)N^l\  (2.8) 

and  whose  t  components  have  subscripts 

m  =  mt_1(l)Aff|\  . mi(f).AfW. 

(Recall  that  the  subsequence  <mi>  depended  on  the  level  L.)  This  completes  our 
recursive  definition. 

Now  we  come  to  the  most  important  definition  of  this  paper. 

DEFINITION  1.  Let  G{,+1,  l  =  0,  1 . A-l,  be  an  Z  + 1-st  level  grid,  occupying  an 

i-th  level  rectangle  Blv,  whose  mesh  points  are  as  given  above.  Let  f  be  any  time 
such  that 

(2-9) 

and  let  f^+1  be  the  greatest  i  +  l-st  level  time  division  point  not  exceeding  t,  An 
L-rJ-st  level  refinement  at  time  t,  corresponding  to  Blv  or  G{,M ,  is  a  sequence  of 
ordered  pairs 

where  j  has  the  values  (2.B).  The  first  components  comprise  the  sequence  of 
i  +  l-st  level  spatial  division  points  contained  in  the  horizontal  sides  of  the 
refinement  rectangle  B\,  (equivalently,  the  sequence  of  x  components  of  the 
grid  points  in  Gt+1);  the  second  components  are  the  approximate  solution  values 
evaluated  at  these  spatial  points,  but  at  time  t^1.  Here  vj+l(t)  is  an  approxima¬ 
tion  to  the  vector  it  (xj  41 .  t ). 
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An  important  property  of  our  definition  is  that  an  i  +  l-st  level  refinement 
exists  not  only  at  i  +  l-st  level  time  division  points  T1*1,  but  also  at  "finer"  time 

division  points  Tuz . Ts  satisfying  (2.9).  However,  solution  values  for  an 

1  +  i-st  level  refinement  are  only  updated  at  i  +  1-st  level  time  division  points. 

For  l  7>  1  let  Blv  be  any  refinement  rectangle,  and  R\fx  its  corresponding 
refinement.  A  vertical  side  of  Blv  which  does  not  lie  on  the  boundary  of  the 
region  R  will  be  called  a  coarse /fine  interface.  Similarly,  the  left  or  right  end¬ 
point  of  R[fx  will  also  be  called  a  coarse/ftne  interface  if  it  does  not  lie  on  the 
left  or  right  boundary  of  the  region  R. 

The  first  level  (or  coarse)  space-time  grid  occupies  the  whole  rectangle 
Bf  =  R.  Hence,  the  first  level,  or  coarse,  refinement  is  present  at  all  times,  and 
higher  level  refinements  are  considered  to  be  superimposed  on  it.  (Strictly 
speaking,  we  should  not  call  this  a  refinement,  since  it  doesn't  refine  anything. 
We  use  this  terminology  to  avoid  special  cases.)  Wc  will  assume  as  given  the  larg¬ 
est  wave  propagation  speed.  This  is  usually  known  by  the  problem  originator, 
and  determines  the  spacing  of  the  coarse  refinement. 

Another  factor  which  must  determine  the  spacing  of  the  coarsest 
refinement  is  the  wavelength  of  any  "background  disturbances"  to  the 
phenomenon  of  interest  (see  our  model  problem  PI  later  in  this  chapter  for  an 
example).  This  too  is  assumed  known;  for  guides  to  the  number  of  mesh  points 
needed  per  wave  length,  see  Kreiss  and  Oliger  [29]. 

We  will  now  discuss  some  further  restrictions  imposed  on  our  refinement 
rectangles.  We  will  require  that  no  two  i-th  level  refinement  rectangles  in  the 
same  i- th  level  horizontal  strip  can  intersect  or  abut.  (But  i-th  level  rectangles 
in  adjacent  strips  may  abut.)  Assume  an  1-th  level  strip  contains  two  i-th  level 
rectangles  BlM  and  Blv,  having  left  and  right  vertical  sides  with  x  coordinates 


13 


*«(#*)  •  xu[n)  *«(»')  •  Xu{v) ' 

respectively.  Without  loss  of  generality,  assume  that  the  left  side  of  the  former 
is  to  the  left  of  the  left  side  of  the  latter,  a(fi)  <  a(i/).  Then 

ufjx)  <  a(u). 

This  is  no  restriction  in  practice;  if  two  such  rectangles  overlap  or  abut,  we  sim¬ 
ply  consider  them  to  be  one  rectangle. 

Let  \i  ~  kt/  Then  our  construction  ensures  that  X<  is  a  constant  depend¬ 
ing  on  l .  For  simplicity,  our  implementation  restricts  the  refinement  ratios  for 
l  as  l  to  be  the  same.  i.e.,  =  N  and  M ^  -  M,  for  l  =  1,  2 . A-l.  This  con¬ 

dition  is  not  essential,  but  it  poses  no  real  restriction,  as  we  will  see  in  Section 
8.3.  For  convergence  studies,  we  shall  in  addition  assume  that  M  =  N,  so  that 
\t  =  constant,  independent  of  l . 

In  Sections  4  and  5  we  assume  that  we  adjust  the  mesh  only  at  coarse  time 
steps.  Thus  all  partitions  Pt  are  the  same  as  the  first  level  partition.  In  prac¬ 
tice,  we  can  adjust  the  mesh  more  often  (or  less  often)  than  every  coarse  time 
step,  but  for  each  l  we  always  adjust  at  times  which  are  a  multiple  of  the  1-th 
level  time  step.  That  is,  every  partition  Pt  is  of  the  form 

0  =  ^  <  •  •  •  <  “  T, 

where  —  i5(£)  is  a  small  positive  integer. 

The  blackened  rectangle  in  Fig.  1  illustrates  the  possibilities.  If  we  adjust 
the  mesh  only  at  coarse  time  steps,  then  this  rectangle  contains  no  mesh  points 
in  its  interior.  However,  if  we  allow  mesh  adjustment  between  coarse  time  steps, 
then  the  blackened  rectangle  may  contain  six  subrectangles 

2.3.  Operations  on  Refinements.  We  now  describe  the  operations  we  can 
perform  on  refinements.  For  simplicity,  we  shall  assume  that  we  adjust  the 
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mesh  only  at  coarse  time  steps,  and  use  the  notation  (2.6)  (or  the  division  points 
between  strips. 

We  shall  say  that  two  refinements  are  equivalent  when  their  first  com* 
ponents  (x  coordinates)  are  the  same,  regardless  of  the  time  or  the  solution 
values.  Thus,  for  all  times  (2.9)  encompassed  by  the  refinement  rectangle  Blv, 
the  refinements  Rlv*l(t)  corresponding  to  Blv  are  equivalent.  In  this  sense,  we 
can  say  that  to  each  rectangle  Blv  or  grid  Gi+1  there  corresponds  one 
refinement.  This  equivalence  concept  is  useful  for  describing  refinement  mani¬ 
pulations  which  do  not  depend  on  the  differential  equation  calculations.  Clearly, 
only  refinements  with  the  same  level  can  be  equivalent. 

Suppose  first  that  there  is  an  i-th  level  refinement  rectangle  B^  {l  >  0)  con¬ 
tained  in  the  strip  (2.6).  Assume  that  the  horizontal  sides  of  B1^  occupy  the 
interval 

xa{/i)  Z  £  xu{p,)- 

Also  assume  that  no  part  of  any  i-th  level  refinement  rectangle  in  strip  St_i  lies 
in  this  interval.  Then  we  will  say  that  the  refinement  RlJl  corresponding  to  B\ 
has  been  created  at  time  t  =  t'~l.  Similarly,  if  we  replace  5t-i  by  5t+i  and  ft_1 
by  i\  we  say  that  Rl*x  has  been  deleted  at  time  f*. 

Now  suppose  there  are  two  i-th  level  refinement  rectangles  Bj^cSi  and 
BlcSi+t.  According  to  our  definition,  the  refinement  RlJl  corresponding  to  B J* 
only  exists  for  f*'1  s  i  if*,  and  the  refinement  R\?x  corresponding  to  Blv  only 
exists  for  f*  £  t  £  fiM.  We  will  now  examine  the  possible  relationships  between 
these  refinements. 

Suppose  the  rectangles  B\  and  Blv  have  the  same  left  and  right  sides, 
a{fu)  =  a(u)  and  uQu)  =  u(v).  Then  the  first  components  of  the  refinements 
corresponding  to  these  rectangles  are  the  same.  By  our  definition,  the 
refinements  corresponding  to  B and  B\  are  equivalent.  In  this  sense  we  may 


say  that  a  single  refinement  now  exists  for  times  Zt-1  ti*1. 


Now  suppose  that  the  refinement  rectangles  are  situated  as  before,  but 

a[fj)  *  a(t/)  <  «(i/)  < 

(with  at  most  one  equality),  and  no  part  of  any  other  Z-th  level  refinement  rec¬ 
tangle  in  strip  Si+1  lies  in  the  interval 

*«(m)  *  *  *  *  Loo- 

Then  we  will  say  that  the  refinement  Rlpl  has  contracted  at  f  =  f ‘  to  form  the 
refinement  R1*1 .  By  interchanging  refinement  rectangles  and  strips,  respec¬ 
tively.  an  analogous  definition  can  be  given  for  an  expanding  refinement. 

If  and  Blv  are  situated  as  before,  but 

a(fi)  <  a(i/)  <  u(jx)  <  w( v). 

and  no  part  of  any  other  Z-th  level  refinement  rectangle  in  strips  St  or  St+, 
occupies  the  interval 

*Lo*)  . 

then  refinement  Rjf1  has  moved  right  at  f  =  f  *  to  become  the  refinement  R\fx . 
Analogously,  we  can  define  what  it  means  for  a  refinement  to  move  left. 

Finally,  suppose  rectangle  Bp  is  in  strip  St  as  before,  but  strip  St4. j  contains 
two  (disjoint)  Z-th  level  refinement  rectangles  Blv  and  Bln,  with  the  former  to  the 
left  of  the  latter.  Assume  that 

a(v)  s  a{fx)  <  u(v)  <  a(rr)  <  o(fi)  s  «(tt), 

and  that  no  part  of  any  other  Z-th  level  refinement  rectangle  in  strips  or  Si+t 
lies  in  the  interval 


*L(v)  *  *  * 

Then  the  refinement  R1*1  is  said  to  separate  or  split  into  refinements  /?*,*’  and 


R**1 .  Analogous  definitions  cam  be  given  for  two  refinements  to  merge  into  a 
third. 

The  above  are  typical  operations  on  refinements,  but  they  do  not  exhaust 
the  possibilities  (for  example,  a  refinement  could  split  into  three  refinements, 
although  this  is  quite  rare).  Fortunately,  however,  an  exhaustive  listing  is  not 
needed.  All 'that  is  required  is  aui  algorithm  which  takes  a  set  of  i-th  level 

refinements  (i  >  1)  at  time  t  =  ti,  i  =  0,  1 . s-1  and  produces  a  new  set  of 

such  refinements.  For  each  i,  once  the  left  and  right  edges  of  the  new 
refinements  are  determined  (by  local  error  estimates),  this  readjustment  can 
be  done  in  a  single  left-to-right  scan  of  the  existing  1-th  level  refinements. 

In  Section  6  we  will  explain  how  these  operations  are  implemented. 

2.4.  A  Model  Problem.  We  will  now  introduce  our  first  model  problem,  which 
will  be  used  in  some  of  our  computations  in  Section  0.  It  is  the  first  order  wave 
equation  ("color  equation") 

Ut  =  -cu*,  aszsb,0st,0<cl  (Pi) 

u(x,0)  =  0(x),  a  £  z  £  b , 

u(0,t)  =  g(-ct),  0  js  f , 

with  exact  solution  u{x,t)  -  p(x-cf).  We  take  a  =  0,  b  =  4,  and  c  =  1.  The 
function  g  is  taken  to  be  a  Gaussian  pulse,  traveling  to  the  right  with  speed  c, 
superimposed  on  a  sinusoidal  background, 

g{x)  =  exp(-a(z+^)z)  +  0.1sin2fr(x+fc), 

with  a  =  200.  The  parameter  a  control  the  steepness  and  thickness  of  the  pulse. 
For  a  =  200,  the  pulse  occupies  about  0  percent  of  the  interval  [0,  4].  This 
models  more  realistic  problems  such  as  an  atmospheric  front  or  storm  surge. 


s.\- 
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2.5.  Solution  Algorithm.  We  now  describe  our  algorithm.  As  in  the  initial 
value  problem  for  ordinary  differential  equations,  we  will  need  to  give  a  toler¬ 
ance  6  on  the  local  truncation  error,  which  will  be  used  to  decide  where  to  refine 
the  mesh.  (We  have  only  used  absolute  error  since  all  of  our  example  problems 
vary  between  0  and  1.  In  general,  one  should  use  a  combination  of  relative  and 
absolute  error,  as  in  Shampine  and  Gordon  [42].) 

For  the  initial  value  problem  for  o.d.e.’s,  Stetter  [46],  and  others,  have 
shown  how  to  estimate  (but  not  control)  the  global  error  while  the  solution  is 
being  computed.  This  requires  only  a  small  amount  of  additional  computation 
and  storage  for  that  case.  Further  investigation  would  be  needed  to  apply  this 
to  the  initial  boundary  value  problem.  But  even  if  were  done,  one  would  still 
need  to  prescribe  a  local  error  tolerance. 

Since  our  refinements  are  defined  recursively,  we  can  define  the  algorithm 
recursively.  Wc  will  use  an  Algol-style  notation.  A  non-recursivc  description  of 
the  same  algorithm  on  a  model  problem  was  given  in  [5]. 


procedure  advance. solution(i,  f ); 
real  t ;  integer  1,  i; 
value  t ; 

comment  advance  solution  from  t  to  t  +  M^~^kt\ 
for  i  =  1  to  do 

begin 

if  mesh  adjustment  time  at  this  level 
begin 

estimate. error(l); 

if  1  <  max. level  then  adjust.  mesh(i  +  1) 
end; 

if  level  1  +  1  exists  then  advance. solution^  +  1,  t), 
compute  the  solution  at  time  t  +  kt  on  all  level  1  refinements; 
Set  t  «- 1  +  fct 
end; 

if  l  >  1  copy  (project)  solution  values  from  this 
refinement  to  its  parent  refinement  at  level  1-1; 
end  advance. solution; 

1  =  1; 

1  =  0. 

Obtain  initial  values  on  coarsest  mesh  (level  1); 


linn,  JmnA 


advance . solution^ .  t), 

Thus,  we  advance  on  the  highest  level  refinements  first,  then  the  next 
highest,  etc.  (One  can  also  proceed  from  the  coarsest  level  to  the  finest,  and 
this  may  be  advantageous  in  more  space  dimensions.)  For  example,  if  our  prob¬ 
lem  has  three  refinements,  (one  each  at  levels  1,  2,  and  3,  respectively)  and 
refinement  ratios  N  =  M  =  3,  we  compute  the  solution  at  the  next  coarse  time 
step  by  advancing  the  different  levels  in  the  following  order:  3,  3,  3,  2.  3.  3.  3,  2. 
3,  3,  3,  2,  1.  (Notice  that  each  recursive  activation  of  our  procedure  has  its  own 
private  copy  of  i  and  t,  which  may  be  different  from  the  values  of  i  and  t  at 
other  levels.  Thus  t  must  be  called  by  value  and  not  by  name  or  reference.) 

We  now  explain  the  steps  of  the  algorithm  in  more  detail. 

"Mesh  adjustment  time"  means  that  the  time  level  is  one  of  those  occurring 
in  the  partitions  Pt  given  in  Section  2.2.  If  the  time  t  is  a  member  of  one  of  the 
partitions  Pt,  we  estimate  the  local  truncation  error  (see  Section  3)  that  would 
be  made  if  we  took  one  forward  time  step  in  the  level  l  refinement,  but  we  do 
not  actually  take  the  step.  Mesh  points  whose  error  estimate  exceeds  the  local 
error  tolerance  are  marked  as  needing  refinement.  These  points  are  grouped 
into  intervals.  Several  extra  "buffer"  mesh  points  are  added  to  both  ends  of 
each  such  interval.  This  will  be  explained  later. 

In  general,  there  may  be  more  than  one  refinement  at  each  level  (except 
the  first).  In  that  case,  the  operations  are  done  for  all  refinements  on  a  given 
level,  starting  with  the  leftmost  refinement. 

To  adjust  the  mesh,  we  compare  the  intervals  just  produced  with  the  exist¬ 
ing  refinements.  If  these  are  not  identical,  refinements  may  have  to  be  "moved", 
created,  deleted,  merged  or  separated.  If  a  refinement  on  level  l  moves  into  a 
region  formerly  occupied  only  by  a  refinement  on  level  l— l,  we  may  need  solu¬ 
tion  values  that  do  not  yet  exist  at  mesh  points  in  refinement  l  These  are 
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obtained  by  linear  or  quadratic  interpolation  in  space  from  solution  values  on 
the  i-l-st  level  refinement. 

Creation  of  a  new  refinement  is  done  the  same  way,  by  spatial  interpolation 
from  its  parent  refinement.  At  any  mesh  adjustment  time,  an  i-l-st  level 
parent  refinement  can  give  birth  to  any  number  of  i-th  level  refinements,  but  no 
higher  level  ones.  An  exception  is  made  at  t  -  0.  If  refinement(s)  of  the  coarse 
mesh  are  needed  at  that  time,  we  obtain  the  new  solution  values  directly  from 
the  initial  function  /  rather  than  from  interpolation.  That  is.  we  take  as  initial 
values 

*J(0)  */(«+#»).  (210) 

where  j  =  0,  1 . JV(0)  when  i  =  1,  and  j  has  the  values  given  in  (2.B)  for  any 

other  refinement  introduced  at  t  =  0.  This  allows  us  to  add  as  many  levels  of 
refinement  as  are  initially  necessary.  Thus,  the  method  performs  properly  even 
when  the  initial  mesh  is  "too  coarse". 

If  a  refinement  occupies  a  spatial  interval  /,  it  can  be  deleted  when  it  has 
no  children,  and  the  local  error  estimate  of  its  parent  in  interval  /  is  below  the 
tolerance. 

We  now  explain  the  "copy"  operation.  Certain  points  (*,  0  of  our  region  R 
(such  as  F  in  Fig.  2)  are  covered  by  mesh  points  in  refinements  at  more  than  one 
level.  (This  is  done  for  ease  of  implementation,  and  arises  because  we  overlay, 
rather  than  patch  or  interpolate  refinements.)  We  always  want  to  use  the  solu¬ 
tion  values  computed  on  the  finest  (highest  level)  mesh,  so  for  these  points  we 
must  project  (copy)  the  values  obtained  on  a  level  l  refinement  to  the  appropri¬ 
ate  positions  in  the  parent  (level  i- 1)  refinement. 

It  remains  only  to  describe  our  difference  approximations. 
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2.8.  Difference  Approximations.  We  first  define  the  forward,  backward,  and 
centered  difference  operators  Di ,  and  ,  operating  on  i-th  level 
refinements: 

D'+vJLt)  =  hc\E  -  1  )vv(t)  =  -  vv(t ))/ht. 

Divv(t)  =  V‘(l  “  E~l)vv(t)  =  (vv(t)  -Vy-iitfi/hi, 

DM)  =  (2Aj)-‘(^  -E~')vv(t)  =  («„*,(0  -*„_,(*))/ 2A*. 

where  we  have  omitted  the  superscript  Z  on  v. 

To  advance  the  solution  one  time  step  from  t  to  f  +  ifct,  we  will  treat 
separately  the  interior  of  a  refinement,  coarse /fine  interfaces,  and  boundaries. 
In  the  interior  of  an  i-th  level  refinement  we  will  use  explicit  two  (time)-level 
difference  approximations  to  (2.1), 

vl(t+kt)  =  Q0vlM)  +  *»n(0.  (2-11) 

where  t  -  tlm. 

<?o=<?o(0=  t  *.*)£*. 

J*-r 

E  =  E(l)  is  the  shift  operator 

fy/t(0  s  vUj(0. 

g  and  r  are  nonnegative  integers,  vlv(t)  is  an  approximation  to  u(xlv,  t ),  and 
/t(<)  =  F{xlv,  t).  (By  the  interior  of  a  refinement,  we  mean  all  its  points  except 
the  r  leftmost  and  q  rightmost  ones.)  The  coefficients  Ai  are  assumed  to  depend 
smoothly  on  their  arguments.  For  example,  the  interior  Lax-Wendroff  approxi¬ 
mation  to  our  model  problem  (with  t  -  tlm  =  ™Jct)  is 

vj{t+kt)  =  (/  -  cktDl0  +  fa%zDl,Di)vj(t). 

The  restriction  to  two-level  schemes  is  necessary  to  simplify  manipulations 
with  refinements.  (When  the  spatial  mesh  is  adjusted  at  time 
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t*,i  -  0.  1 . *— 1,  it  would  be  awkward,  and  require  more  storage,  to  adjust 

the  mesh  at  previous  time  levels  too.)  This  also  simplifies  error  estimation,  as 
will  be  seen  in  Section  3.  (This  restriction  does  not  exclude  two-level  schemes 
with  fractional  time  steps,  such  as  two-step  Lax-Wendroff.) 

The  restriction  to  explicit  schemes  is  more  fundamental.  This  is  no  restric¬ 
tion  for  purely  hyperbolic  problems,  but  can  be  a  restriction  for  more  compli¬ 
cated  problems  ( e.g . ,  coupled  heat  and  sound).  Since  our  algorithm  calculates 
solutions  at  a  given  time  level  piecewise  in  various  parts  of  the  interval 
asis6,  this  restriction  is  not  merely  for  convenience. 

In  Section  3  we  will  also  require  that  the  local  truncation  error  (per  unit 
time  step)  of  the  interior  approximation  must  have  the  same  order  in  space  and 
time.  Since  we  will  most  often  use  interior  approximations  which  are  second 
order  in  space  and  time,  this  restriction  is  not  too  severe. 

At  coarse-fine  interfaces  (between  level  l  +  l  and  level  l  refinements)  which 
do  not  abut  boundaries,  we  use  a  hybrid  method,  the  coarse /fine  approximation 
[9]  corresponding  to  difference  method  (2.11).  We  apply  the  right  side  of  (2.11) 
to  solution  values  on  level  l  with  space  step  Aj  and  time  step  ifc{+ j,  for 
i  =  1,  2,  ■  •  •  .  M  (that  is,  replace  \  by  tfcj+i/hJ,  and  store  the  solution  value  on 
the  level  l  + 1  refinement. 

For  Lax-Wendroff  applied  to  our  model  problem  this  is 

■****♦!)  =  V  +  J £*ic*&|£U7!.)vj(im)- 

Fig.  2  shows  the  stencil  at  the  interface  between  a  refinements  on  levels  2  and  3 
(f  =  2)  when  N  -  M  -  3,  so  that  the  third  level  space  step  /ts  =  hz/  3  and  the 
third  level  time  step  fc3  =  *2/3.  When  i-  1,  points  A,  B,  and  C  are  used  to 
advance  to  point  D.  When  i  =  2,  points  A,  B,  and  C  are  again  used  to  advance  to 
point  E,  and  when  i  =  3,  the  same  points  are  used  to  advance  to  point  F. 


v“  vS  O  V 
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For  a  difference  scheme  (2.11)  whose  stencil  is  more  than  three  mesh 
points  wide,  we  will  need  to  use  the  coarse/flne  approximation  r  times  in  Fig.  2 
(and  correspondingly  g  times  at  the  right  end  of  a  refinement).  This  is  done. 
e.g. ,  for  r  =  2,  by  using  the  stencil  illustrated  in  Fig.  2  to  get  point  D.  then  shift¬ 
ing  the  stencil  one  finer  mesh  point  to  the  right  to  get  the  point  to  the  right  of  D. 
(This  involves  a  spatial  interpolation  in  the  coarser  mesh,  which  is  done  as  in  the 
mesh-adjustment  step). 

Finally,  boundaries  are  treated  the  same  as  with  a  uniform  mesh;  if  an  l- th 
level  refinement  abuts  the  left  boundary,  we  set 

vlM(t+k,)  =  £  S^v^t-aki)  +gl(t).  M  =  0.1 . r- 1,  (2.12) 

where 

s^(l)  =  £  CfrHzl+jhi.  t-oki.ht)#,  a  =  0.-1 

t  -  tjn  s  miei  ,r<r,  and  =  0.  The  approximation  at  the  right  boundary 

is  analogous.  For  our  model  problem  we  use  the  prescribed  values 

vj(t+kt)  -  g  (-c  (t  +fcj)) 

at  the  left  boundary;  and  upwind  differencing 

vi(f+*j)  -  (I  ~  cktDl-)vj(t) 

at  the  right  boundary. 

Once  again  we  have  restricted  ourselves  to  two  time  levels,  for  the  same 
reasons  as  before.  We  allow  ourselves  implicit  boundary  conditions  here 
(S^  *  0)  since  we  can  first  solve  for  the  points  on  the  right  hand  side  of  (2.12) 
using  the  explicit  interior  approximation. 

An  extremely  important  feature  of  this  method  is  the  use  of  a  buffer  on 
either  end  of  any  refinement  (except  the  coarsest  one).  If  we  are  estimating  the 
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truncation  error  for  the  refinement  j xjf  and  the  error  tolerance  is  exceeded 
between  j  =  a  and  j  =  u,  then  we  instead  refine  from  j  =  a  -  6t  to  j  -  a  +  bt , 
where  b4  is  the  buffer  length  for  refinements  of  level  Z+l.  That  is.  both  ends  of 
the  Z  +  l-st  level  refinement  are  padded  with  bt  extra  cells  of  width  V  In  gen¬ 
eral.  if  our  1-th  level  refinement  requires  several  intervals  of  1  +  1-st  level 
refinement  (according  to  the  error  estimate),  then  each  such  interval  is  padded 
as  above.  (This  may  cause  some  1  +  1-st  level  refinements  to  merge.) 

How  do  we  choose  b{?  From  Fig.  2,  we  see  that  bt  should  be  at  least  one. 
because  we  use  the  coarse /fine  approximation.  For  safety  we  make  it  two.  We 
shall  assume  that  the  (given)  maximum  wave  speed  is  c .  For  simplicity  we  shall 
assume  that  all  partitions  fl.la  1.  are  the  same,  and  that  we  check  the  error 
every  coarse  time  steps.  Therefore,  in  time  kx  a  wave  could  travel  left  or  right 
a  distance  of  d}kx  =  ciSXiAi  =  ci5 \xNl~lhi,  or  ci>X1JVt-1  cells  of  width  hj.  (Here 
the  Z-l  is  an  exponent,  not  a  superscript.)  So  we  take  b{  ~  2  +  IciJX,lVl-1V  where 
fxl  is  the  ceiling  function  (the  least  integer  greater  than  or  equal  to  x).  (For 
difference  approximation  (2.11),  bt  must  be  modified  by  replacing  2  by  g+1  at 
the  left  end  of  a  refinement,  and  by  r+1  at  the  right  end.)  Obviously,  higher  level 
refinements  have  larger  buffers. 

The  buffer  mechanism  has  several  beneficial  consequences.  First,  and  most 
important,  it  insures  that  the  rapidly  varying  part  of  the  solution  does  not 
escape  into  the  coarser  region.  As  we  saw  in  Section  1,  this  is  absolutely  essen¬ 
tial  to  the  success  of  the  algorithm.  Secondly,  this  policy  allows  us  to  use 
difference  approximations  at  coarse /fine  interfaces  which  would  otherwise  not 
be  accurate  enough  in  the  fine  mesh.  We  can  also  estimate  the  local  truncation 
error  at  coarse /fine  interfaces  in  a  very  simple  manner  (see  Section  3.2).  Third, 
it  allows  "smooth"  transitions  in  mesh  width.  That  is,  a  level  l  refinement  can 
abut  a  level  l  +  l  or  l-l  refinement,  but  not  others.  This  is  important  when  using 
recursive  refinements.  Fourth,  it  keeps  the  refinements  from  splitting  into  tiny 
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pieces,  because  level  i  +  1  level  refinements  which  are  closer  than  26(  level  l  cells 
apart  (before  buffering)  are  joined  together.  (If  the  local  truncation  error  were 
large  in  absolute  value  but  suddenly  changed  sign,  this  might  cause  splitting  into 
pieces.)  We  make  this  condition  even  more  stringent  by  joining  together  any 
level  l  +  l  refinements  which  are  less  than  2bt+2  level  l  cells  (of  length  ht)  apart 
(before  buffering).  Fifth,  buffering  allows  us  to  specify  a  priori  the  times  to 
check  the  local  error  (and  adjust  the  mesh).  In  particular,  we  need  not  check 
the  error  at  every  (fine)  time  step  (Section  8  shows  that  this  is  very  expensive). 
We  can  instead  check  at  every  coarse  time  step,  or  even  every  coarse  time 
steps,  where  i?  is  a  small  positive  integer.  Sixth,  buffering  contributes  greatly  to 
the  robustness  of  the  algorithm.  Buffers  make  the  algorithm  relatively  insensi¬ 
tive  to  small  inaccuracies  in  the  local  error  estimation. 

Let  us  comment  on  the  storage  required  for  this  algorithm.  If  we  use  a 
two-level  method  and  perform  the  operations  in  the  order  given,  then  we  need 
two  levels  of  solution  values,  just  as  for  a  uniform  mesh.  As  soon  as  all  the  solu¬ 
tion  values  in  a  refinement  at  a  new  time  level  are  known,  we  can  overwrite  them 
on  the  old  solution  values.  A  slight  amount  of  additional  storage  is  needed  for 
pointers  and  indices;  this  is  minuscule  compared  to  space  for  solution  values. 
Next,  free  space  is  needed  to  separate  the  solution  values  on  refinements.  The 
amount  is  variable,  but  can  be  chosen  quite  small  (see  Section  6.3).  Finally, 
storage  is  needed  for  error  estimates;  but  these  can  be  done  a  refinement  at  a 
time,  so  we  only  need  two  vectors  (when  solving  a  scalar  equation),  each  the  size 
of  the  largest  refinement.  We  did  not  implement  our  algorithm  in  a  way  which 
minimizes  the  amount  of  storage. 

3.  Estimation  of  the  Local  Truncation  Error.  In  this  section  we  show  how  to 
estimate  the  local  truncation  error.  We  examine  two  methods  for  estimation  in 
the  interior  of  refinements:  differences  and  Richardson  extrapolation.  Then  we 
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propose  methods  for  estimating  the  error  at  coarse /fine  interfaces  and  at  boun¬ 
daries. 

The  first  important  conclusion  of  this  section  (and  of  the  computations  in 
Section  6)  is  that  several  methods  can  be  used  to  estimate  the  error,  both  in  the 
interior  and  at  boundaries.  (In  [5]  we  examined  two  other  methods,  neither  of 
which  is  as  useful.)  Thus  our  algorithm  is  quite  general  as  well  as  robust.  The 
second  important  finding  is  that  the  Richardson  method  is  more  convenient  in 
the  interior,  and  differences  are  more  convenient  at  boundaries. 

For  simplicity,  we  shall  write  all  approximations  as  occurring  on  a  uniform 
mesh.  We  will  always  assume  that  \  =  k/h  =  constant,  and  that  the  solution  of 
the  differential  equation  has  sufficiently  many  derivatives.  When  we  speak  of 
asymptotic  estimates  and  leading  terms,  we  shall  always  mean  "as  h  -*  0". 

3.1.  Two  Methods.  We  now  explain  two  methods  for  estimating  the  leading 
terms  of  the  asymptotic  expansion  of  the  local  truncation  error.  These  methods 
require  extra  computations,  in  addition  to  those  required  to  compute  the  solu¬ 
tion.  We  first  examine  the  interior  of  a  refinement. 

3.1.1.  Differences.  We  will  illustrate  this  method  on  an  example.  The  local 
truncation  error  for  the  Lax-Wendroff  approximation  to  problem  PI  on  a  uniform 
mesh  with  stencil  centered  at  (x,f )  is 

^{-khittt{x,t)-hh^{x,t))  +  0(/i4).  (3.1) 

We  use  the  differential  equation  to  replace  t  derivatives  by  x  derivatives  in  this 
expression.  We  then  obtain 

ck ^»(x,f)(cz\2  -  1)  +  0(h*). 

We  may  now  approximate  it m  by  a  five-point  divided  difference  at  points  in 
the  interior  of  a  refinement.  Specifically,  with  t  dependence  omitted. 


XLj (x)  »  (- Zu{x-Zh )  +  Au(x-h)  -4it(x+/i)  +  2u(x+2h))/4/i3. 


3.1.2.  Richardson  Extrapolation.  We  shall  apply  this  method  (which  was 
suggested  by  J.  Oliger)  to  our  linear  hyperbolic  operator  (2.1).  In  the  interior  of 
a  refinement  we  approximate  this  system  by  any  linear  multi-(time)  level  expli¬ 
cit  difference  scheme  whose  local  truncation  error  per  unit  time  step  has  the 
same  order  p  in  space  and  time: 

vv{t+k)  =  Q(h)vv(t)  s  £  Q„vv{t-ak)  +  kFv{t),  v  =  r,  r+1 . N0-q.  (3.2) 

a=0 

Here  x=xlJ  =  a  +  uh,t=tm=  mJc , 

Qo=  t,  Aj'(xv+jh,  t-ak,  h.)E*  a  -  0.  1 . p.  (3.3) 

i-~r 

the  Aj0  are  matrix  coefficients,  and  E  is  the  spatial  shift  operator. 

The  local  truncation  error  of  method  (3.2)  applied  to  (2.1)  at  the  point 
(z.  0  =  (*w.  tm)  is  given  by 


u{x,t+k)  =  £  Qau(x,t  -ok)  +  kF{x.t)  +  k{h*Tx{x,t)  +  k*  [7,(x ,f )) 

(7=0 

+  0{h?+z), 


(3.4a) 


where  Tx  and  Ux  are  sufficiently  smooth  functions  of  z  and  t .  Symbolically, 

u(x,t+k)  =  9(h)u(x,f)  +  kF{x,t)  +  kr(x,t).  >3. 4b) 

For  the  global  truncation  error  e(x,t )  =  u(x,f)  -  vv{t)  we  obtain  the  expansion 

e{x,t+k)  =  9(/i)e(z,f)  +  k  r(z ,  t )  (3.5) 


by  subtracting  (3.2)  from  (3.4).  We  shall  assume  that  the  difference  approxima¬ 
tion  is  defined  for  ail  (x,f),  not  only  at  one  set  of  mesh  points. 

To  use  the  Richardson  method,  we  take  one  (time)  step  of  the  approxima¬ 
tion  (3.2)  with  z  =  x„,  t  s  tm,  using  mesh  spacing  h  in  space  and  k  in  time.  (See 
Fig  3  for  a  two-level  scheme  with  q  =  r  =  1,  with  stencils  shifted  horizontally  for 
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clarity.)  We  then  repeat  the  step  using  (3.2)  with  tm  replaced  by  tm+i  and  the 
same  mesh  spacing,  obtaining  the  approximation  vv(t+2k).  (Before  performing 
this  second  step,  we  will  need  to  generate  more  points  on  level  t  =  (rn  +  l)fc  by 

applying  (3.2)  with  ( xv,tm )  replaced  by  (x^,tm),  j  -  -r.-r+l . 0,1 . g. 

Thus  in  practice  this  estimation  is  done  for  all  interior  points  of  a  refinement  at 
once.)  Then  we  start  over  at  (x,t)  =  (x„,tm)  using  (3.2).  but  with  stencil  spacing 
2h  and  2k  (i.e. .  replace  jh  by  2jh  and  k  by  2k  in  (3.2)  and  (3.3)).  obtaining  the 
approximation  (f  +2k).  We  then  subtract  the  two  approximations  obtained 
at  level  fm+. 2,  and  divide  by  2*+1  -  2  to  obtain  the  desired  estimate. 

We  shall  prove  the  validity  of  this  method  when  both  the  differential  and 
difference  equations  have  constant  coefficients  (they  may  not  depend  on  x  or  f 
but  those  of  the  difference  approximation  may  depend  on  h  when  B  in  (2.1)  is 
nonzero).  However,  our  numerical  results  in  Section  8  will  show  that  this  pro¬ 
cedure  is  applicable  in  much  more  general  circumstances.  In  practice,  our  inte¬ 
rior  difference  approximation  will  always  be  two-level  (p  =  0).  Here  we  do  not 
rewrite  f  derivatives  in  the  local  truncation  error  in  terms  of  x  derivatives. 

THEOREM  1.  Approximate  the  hyperbolic  operator  (2.1)  for  the  Cauchy  prob¬ 
lem  by  the  consistent  multilevel  explicit  interior  difference  scheme  (3.2). 
Assume  that  both  operators  have  constant  coefficients.  If  the  undifferentiated 
term  Bu  in  (2.1)  is  nonzero,  assume  that  the  coefficients  (3.3)  in  the  difference 
operator  are  smooth  functions  of  h.  Assume  that  the  local  truncation  error  per 
unit  time  step  t  (3.4)  has  the  same  order  p  in  space  and  time,  that  the  solution 
u  of  the  differential  equation  and  the  global  error  e  =  u  -  v  are  sufficiently 
smooth  functions  of  x  and  t,  and  that  \  =  k/h  =  constant.  Then  we  can  esti¬ 
mate  the  ( lowest  terms  of  the  asymptotic  expansion  of  the )  local  truncation 
error  kr{z,t)  at  the  point  (x,t )  =  ( x„,fm )  in  the  interior  of  a  refinement  using 
the  Richardson  method,  and 


kr(x,t)  =  k(hPTi(x.t)  +  kP  U^x.t))  +  O^*2) 

=  ( vv(t+2k )  -  vfp(t+2k))/  (2P  +  1  -  2)  +  0(/i**2). 

where  Tx  and  Ul  are  sufficiently  smooth  functions  of  x  and  t ,  vv(t  +2k)  is  the 
approximation  obtained  by  applying  one  step  of  (3.2)  urith  t  replaced  by  tm+, 
and  mesh  spacing  h  and  k,  and  v^it  +2k)  is  the  approximation  obtained  by 
applying  one  step  of  { 3.2)  with  (x.f  )  =  (xv,tm)  and  mesh  spacing  2h  and  2k. 

Proof.  The  local  truncation  error  e&^x.t  +2*)  =  u{x,t+2k )  -vf?)(t+2k)  of 
our  double  size  step  is,  to  leading  order  terms,  2P+1  times  the  local  truncation 
error  for  a  single  step.  That  is, 

e^(x,f  +2*)  =  Q(2h)e  (x,t)  +  2P+,A:t(x.O  +  0(/i*+2),  (3.6) 

obtained  by  replacing  h  by  2h  and  k  by  2k  in  (3.5).  (From  this  formula  and 
(3.3)  we  can  see  that  using  a  scheme  with  more  them  two  time  levels  will  entail 
storing  many  previous  time  levels  of  the  solution.  This  would  be  highly  impracti¬ 
cal  in  multidimensional  problems) 

We  apply  two  single  steps  of  the  method  to  the  error  e(x.  t).  Thus  we 
replace  t  by  t  +  k  in  (3.5)  and  substitute  (3.5)  into  the  result.  This  yields 

e{x,t  +2k)  =  Qz(h)e  (x  ,t)  +  jfcT(x,f+fc)  +  Q(h)kr(x,t) 

=  Qz(h)e  (x ,t)  +  kr{x,t)  +  (/  +  0{h))kr(x,t )  +  0(hp*z) 

=  <?{h)e(x,t)  +  2  kr{x,t)  +  0(/i*+2). 

Here  we  have  expanded  the  first  t  in  a  Taylor  series  about  (x,f)  and  have 
expanded  Q  using  the  consistency.  (The  consistency  relations  express  the 
coefficients  of  the  power  series  in  Q  in  terms  of  the  coefficients  of  the 
differential  equation,  hence  the  former  are  indeed  bounded.)  We  subtract  this 
equation  from  (3.6)  to  obtain  the  computable  quantity 

vv{t  +Zk)  —  (x  ,t  +2k)  =  e^(x,f  +2*)  -  e  (x,f  +2A:) 

=  [9(2/i  )  -  Q*(h)]e  {x  ,t)  +  (2*  +  I  -  2)kr{x,t)  +  0(/i*t2) 
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This  will  yield  our  result  if  we  can  show  that  the  bracketed  expression  is 
0(hp *8).  Since  Q(h)  =  I  +  0(h)  it  is  clear  that  ^(A)  and  9(2A)  agree  up 
through  and  including  first  order  terms  in  their  formal  power  series  expansions. 
So  we  must  show  that  the  terms  with  coefficients  hz  in  the  bracketed  expression 
are  0(hp*z).  Now  e  is  the  same  order  0(hp)  as  r  for  the  Cauchy  problem,  since 
we  have  assumed  smoothness  of  solutions  and  error.  Hence  the  terms  in  ques¬ 
tion  are  0(hz)  times  derivatives  of  e .  hence  0(hp*z).  This  completes  the  proof. 

We  have  called  this  a  Richardson  extrapolation  method,  but  we  are  using  it 
in  a  non-traditional  way.  Both  our  method  and  the  traditional  approaches  (for 
o.d.e.'s  and  elliptic  p.d.e.’s)  improve  the  accuracy  of  the  approximate  solution 
by  estimating  the  local  truncation  error.  But  we  use  the  estimate  to  decide 
where  to  refine;  the  traditional  approaches  add  the  estimate  to  the  approximate 
solution.  (Doing  the  latter  would  not  be  useful  to  us.  since  our  estimation  is  not 
being  done  at  every  time  step.)  As  a  consequence,  the  traditional  approaches 
improve  the  order  of  accuracy  of  the  basic  difference  scheme;  our  method  does 
not. 

In  both  error  estimation  methods,  the  quantity  we  control  in  the  interior  is 
not  the  local  truncation  error,  but  the  local  truncation  error  per  unit  time  step 

|r(x,f)|.s  <5, 

where  S  is  the  local  error  tolerance,  and  |  -|.  denotes  the  maximum  absolute 
value  of  all  components  of  an  n-vector  at  the  position  (x,  t ).  Thi.  s  because  one 
power  of  h  in  accuracy  is  lost  in  going  from  the  local  truncation  error  of  the 
interior  approximation  to  the  global  error  [21]. 

Our  computations  in  Section  8  6  of  [5]  show  that  for  our  model  problem  Pi. 
either  of  our  methods  produces  approximately  equally  accurate  estimates  of  the 
local  truncation  error  in  the  interior.  And  clearly,  Richardson  estimates  are 
more  expensive  to  compute  than  difference  estimates  However,  the  Richardson 
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method  is  considerably  more  convenient  to  use  in  the  interior  of  refinements. 
The  difference  method  requires  us  to  explicitly  compute  the  local  truncation 
error,  and  then  to  replace  t  derivatives  by  x  derivatives.  Even  with  a  symbol- 
manipulation  program  like  MACSYMA,  this  can  be  exceedingly  cumbersome  for 
realistic  coupled  systems.  The  Richardson  method  makes  possible  an  (almost) 
automated  approach  to  local  error  estimation.  One  need  only  know  the  order  p 
of  the  method  and  the  factor  2P+1-2  used  to  divide  the  difference  vv  -v®  of 
the  two  approximations  at  time  t  +  2k . 

3.2.  Coarse/Fine  Interfaces.  Let  us  now  discuss  the  modifications  needed 
for  coarse /fine  interfaces  which  do  not  abut  boundaries.  For  concreteness, 
assume  that  an  Z-th  level  refinement  Rt  has  a  descendant  Z+l-st  level 
refinement  which  does  not  abut  the  left  or  right  boundaries  x  -  a  or  x  -  b . 
This  introduces  two  coarse /fine  interfaces,  namely,  the  ends  of  (See  Fig.  2 
for  the  left  end  of  Recall  that,  the  last  time  we  estimated  the  error,  we 

added  enough  padding  or  buffering  (see  Section  2.6)  to  both  ends  of  flt+,  to 
ensure  that  waves  could  not  escape  it,  plus  two  extra  level  Z  (spatial)  cells.  This 
guarantees  that  we  will  not  need  to  refine  the  ends  of  refinement  /?t+l  (unless 
they  abut  boundaries)  and  assures  "smooth"  mesh  transitions.  Since  our  local 
truncation  error  estimates  are  used  only  to  decide  where  to  refine,  we  can  safely 
set  our  estimate  at  the  ends  of  tf{+l  to  zero. 

The  next  question  is  the  choice  of  estimator  at  mesh  points  which  are  one 
(spatial)  mesh  point  on  the  "fine"  side  of  a  coarse /fine  interface,  or  one  mesh 
point  away  from  a  boundary.  (This  is  for  the  case  of  a  stencil  with  three  adja¬ 
cent  spatial  points,  i.e.,  g  =  r  =  1  in  (2.11)  or  (3.2).  In  the  case  where  g  or  r  is 
greater  than  one,  similar  considerations  apply  to  the  g  or  r  points  on  the  "fine" 
side  of  the  interface.)  Fig.  3  shows  that  the  Richardson  method  does  not  yield  an 
estimate  here.  We  also  set  this  estimate  to  zero,  for  the  same  reasons  as  before. 


31 


3.3.  Boundaries.  Let  us  consider  local  error  estimation  at  boundaries.  On 
the  left  boundary,  there  are  J  boundary  conditions  specified  for  the  differential 
equations  (2.1)-(2.4).  We  can  approximate  these  in  the  obvious  way  with  no  local 
truncation  error.  We  will  call  these  "exact"  boundary  approximations. 

When  rilin  the  interior  difference  approximation  (2.11)  or  (3.2),  then  we 
need  n -J  "extra"  boundary  conditions  at  the  left  boundary, 

vM(f+fc)  *  £  Sp>vr(t-ok)  +pM(f).  M  =  0.  (3.7) 

where  S ^  is  as  given  in  (2.12),  but  with  the  appropriate  time  level.  If  r  >  1.  we 

also  need  n(r-l)  additional  boundary  conditions  of  type  (3.7)  for  n  =  1 . r-1. 

(Similar  statements  hold  at  the  right  boundary,  with  J  replaced  by  n -J  and  r 
replaced  by  q.  We  will  only  discuss  the  left  boundary;  the  right  is  similar.)  We 
will  first  consider  the  extra  conditions;  at  the  end  of  the  next  subsection  we  shall 
examine  the  "exact"  boundary  approximations. 

The  local  truncation  error  kr  of  (3.7)  is 

u(*M,f+fc)=  £  Sj^u{xT,t-ak)  +  g^t)  +  kr{x,t),  n  =  0.1 . r-1. 

»*-i 

For  a  restricted  class  of  boundary  approximations,  it  is  tempting  to  recycle 
Theorem  1,  using  boundary  approximations  in  place  of  interior  operators. 
Unfortunately,  this  fails  because  the  boundary  operator  was  not  the  only  opera¬ 
tor  used  to  produce  solution  values  at  previous  time  levels.  (For  example,  if  we 
use  the  first  order  upwind  boundary  approximation  and  the  Lax-Wendroff  inte¬ 
rior  approximation  on  our  problem  Pi  (the  first  order  wave  equation),  then  we 
apply  the  former  three  times  and  the  latter  once  to  obtain  a  boundary  esti¬ 
mate.) 

Thus,  the  Richardson  method  is  not  very  useful  at  boundaries  for  several 
reasons.  First,  wc  must  do  a  complete  error  analysis  of  the  Richardson  method 
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with  both  boundary  and  interior  stencils  for  each  individual  problem.  (We  did 
this  for  problems  PI  and  P2;  see  Section  8.7.)  This  cannot  be  done  in  an 
automated  way.  Second,  there  are  relatively  few  boundary  approximations 
which  have  the  same  order  spatial  and  time  error.  Third,  this  method  does  not 
work  for  implicit  approximations. 

So  we  must  use  differences.  In  contrast  to  the  interior  approximation,  it  is 
usually  practicable  to  write  down  the  local  truncation  error  for  the  boundary 
approximation,  and  rewrite  t  derivatives  in  terms  of  x  derivatives.  For  example, 
in  our  first  order  wave  equation,  if  we  use  upwind  differencing  at  the  right  boun¬ 
dary 

vv(t+Jc)  =i i„(t)  -  c\(v„(t) 
the  local  truncation  error  is 

%-(k2Uu  ~  c  XhhLxx)  +  0(/i3); 
replacing  t  derivatives  by  x  derivatives  yields 

fcfc/i(cX-lX*  +  0(/i3). 

We  then  replace  the  term  by  a  three-point  one-sided  difference. 

The  local  truncation  error  for  extrapolation  boundary  conditions 
(/iZ?+V'i/0(f+/fc)  =  0,  forftxedjil 

can  only  be  estimated  by  differences.  As  an  example,  for  j  =  2.  we  estimate  the 
truncation  error 

+  0(A3) 

by  using  a  four-point  one-sided  difference  (since  the  three-point  estimate  yields 
zero). 

In  Section  8.7  we  numerically  compare  different  methods  of  error  estima¬ 
tion  at  boundaries 
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Gustafsson's  [21]  analysis  and  our  computations  in  Section  6.5  show  that 
the  order  of  accuracy  of  the  boundary  approximation  may  be  one  order  lower 
(but  not  less)  than  that  of  the  interior  approximation,  in  order  to  preserve  the 
global  order  of  accuracy,  which  is  then  the  same  as  the  order  of  the  local  error 
per  unit  time  step  for  the  interior  approximation.  This  means  that  when  we  are 
deciding  whether  to  refine  at  the  boundary,  we  should  not  control  the  local  error 
per  unit  time  step,  but  instead  the  local  error, 

|fcT(z,f)|.s  6. 

3.4.  Boundary  Anomalies.  We  will  now  discuss  an  anomaly  which  arises  fre¬ 
quently  in  model  problems,  but  quite  seldom  in  realistic  problems.  This  prob¬ 
lem  can  only  arise  with  (a)  a  single  scalar  equation;  (b)  annxn  system  which 
can  be  decoupled  into  independent  subsystems;  or  (c)  a  system  with  the  pro¬ 
perty  that  all  the  characteristic  variables  at  a  boundary  are  inflow  variables,  i.e., 
all  the  characteristics  point  into  the  region  (e.g.,  supersonic  inflow).  If  this  is 
the  left  boundary,  then  for  problem  (2.1)-(2.4)  in  diagonal  form,  the  boundary 
condition  (2.3)  has  5n  -  0  (since  un  has  no  components,  ie.,  u  -  u1). 

This  problem  arises  because  no  coupling  occurs  between  inflow  and  outflow 
variables.  Using  the  obvious  difference  approximation  to  these  inflow  boundary 
conditions  in  the  differential  equation  (the  so-called  "exact"  boundary  approxi¬ 
mations)  yields  zero  local  truncation  error. 

Examining  our  problem  Pi,  we  can  see  a  possible  difficulty  with  using  zero 
in  our  error  estimation.  The  left  boundary  condition  contains  a  forcing  (inhomo¬ 
geneous)  term  g  which  results  in  the  wave  entering  the  region.  Clearly,  we  want 
to  put  refinements  around  any  "large"  wave  entering  the  left  as  soon  as  possible. 
If  we  set  the  truncation  error  estimate  to  zero  at  the  boundary,  we  will  not 
detect  the  wave  until  it  has  already  entered  the  region.  This  can  be  remedied  by 
treating  the  forcing  term  g  j  or  gz  of  (2.3)  or  (2.4)  as  generating  a  fictitious  local 
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truncation  error.  For  our  first  order  wave  equation,  where  it(0,i)  =  g{t),  we 
write  down  the  local  truncation  error  (3.1)  for  the  interior  approximation,  then 
replace  x  derivatives  by  t  derivatives,  using  the  differential  equation.  Since 
UtuiO.t)  =  gttt,  we  can  analytically  differentiate  g  to  arrive  at  the  result.  This 
procedure  was  used  in  our  computations  with  problem  Pi  (see  Section  8.4). 
(Notice  that  we  used  an  interior  error  for  a  boundary  approximation,  and  then 
controlled  the  local  error  per  unit  time  step.  We  could  equally  well  have  used 
the  boundary  error,  which  depends  on  utt  and  controlled  the  local  error.  In 
either  case  we  would  control  the  same  power  of  h.) 

In  case  (c),  we  may  also  proceed  as  in  our  first  order  wave  equation.  For  a 
local  truncation  error  of  the  form 

au*z*  +  0^  (3.8) 

we  rewrite  x  derivatives  in  terms  of  t  derivatives  using  the  differential  equation, 

it,  =  A~l{ut  -  Bu  -/), 

and  replace  the  first  term  in  (3.8).  Since  u}u  =  (g  j)iu .  we  can  again  analytically 
differentiate  to  arrive  at  the  result.  In  practice,  the  rewriting  may  be  cumber¬ 
some. 

In  the  vast  majority  of  cases,  our  problem  will  not  have  property  (a),  (b)  or 
(c).  Then  some  of  the  components  of  the  difference  approximation  have  nonzero 
local  truncation  error,  and  our  usual  error  estimates  for  these  component(s) 
will  detect  any  incoming  wave,  since  the  boundary  conditions  are  coupled.  This 
technique  was  used  in  our  problem  P2  (the  second  order  wave  equation)  in  Sec¬ 
tion  8. 

4.  Stability.  In  this  section  we  give  a  brief  discussion  of  stability. 

We  begin  with  a  definition.  The  square  of  the  (discrete)  Iz(x)  norm,  of  an  n- 
vector  v  defined  on  our  refined  grid  at  time  f‘.  t  =  0,  1,  •  ■  ■  ,  s  is  a  sum  of 
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terms,  one  for  each  mesh  point  existing  at  the  given  time.  (If  a  point  (*,  t)  is 
covered  by  more  than  one  grid  point  at  different  levels,  use  the  finest  one.)  Each 
term  is  of  the  form  A{fyj(f<)|f ,  where  is  the  interval  to  the  left  of  the  grid 
point,  and  |  -|2  denotes  the  sum  of  squares  of  the  components  of  the  n-vector. 
(At  the  left  boundary,  use  the  interval  to  the  right  of  the  mesh  point.) 

The  usual  definition  of  stability  for  the  initial  boundary  value  problem  is 
Definition  3.3  of  Gustafsson,  Kreiss  and  Sundstrom  [22]  (hereafter  referred  to  els 
the  GKS  definition).  But  there  are  several  problems  with  this  approach,  if  we 
wish  to  use  this  definition  to  prove  convergence. 

The  GKS  definition  assumes  that  the  initial  data  function  /  is  zero.  This 
may  be  acceptable  for  t  -  0,  but  after  we  integrate  over  the  horizontal  strip  S1( 
we  in  effect  start  a  new  initial  boundary  value  problem  at  f  =  fl,  and  now  the 
"initial"  data  are  not  zero.  It  is  difficult  to  incorporate  a  nonzero  /  into  the  GKS 
definition,  since  it  was  derived  using  Laplace  transform.  However,  it  is  neces¬ 
sary  if  we  wish  to  use  it  to  prove  convergence.  For.  in  an  interval  0  s£  f  £  T  the 
number  of  strips  5(  becomes  unbounded  as  h  -*  0.  The  solution  at  f*  =  T 
depends  on  the  values  of  the  solution  ("initial  data")  at  all  previous  strips,  but 

this  dependence  on  values  at  times  t*~x,  to  t*~z . tx,  has  to  be  removed  if  we 

want  to  prove  convergence.  This  can  only  be  done  if  /  appears  explicitly. 

A  second  difficulty  is  related  to  the  first.  The  GKS  definition  assures  us  that 
exp(-af)  times  the  solution  is  in  iz(x,f)  (the  usual  discrete  Hilbert  space  in  * 
and  t ).  but  for  any  fixed  t  does  not  assure  us  that  the  solution  is  unconditionally 
in  iz(x)  (the  discrete  Hilbert  space  in  x).  (Compare  Theorem  3.1,  GKS.)  In  order 
to  integrate  over  a  new  strip  5<.  we  wish  to  treat  the  solution  values  at  f  =  ft_1 
as  "initial"  values,  and  this  requires  that  they  be  in  fz(x). 

For  these  reasons,  Oliger  [4],  [36]  has  proposed  a  new  stability  definition.  It 
applies  to  approximations  in  any  number  of  space  dimensions,  and  is  the 
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discrete  analog  of  the  well-posedness  condition  for  differential  aquations.  We 
state  it  for  a  problem  in  one  space  dimension. 

DEFINITION  2.  Let  X  =  fej/hj  =  constant,  independent  of  l.  Assume  that  we 
adjust  the  mesh  only  at  coarse  time  points.  The  difference  approximation  (2.10) 
-  (2.12)  on  our  region  R  is  stable  for  a  refined  mesh  (as  described  in  Section 
2.2)  if  for  any  T  >  0,  there  exists  a  constant  >  0  such  that,  for  all  positive 
integers  s ,  all  sets  of  time  division  points  (2.6),  all  kt  >  0. 1  =  1,  2 . A.  satisfy¬ 

ing  our  restrictions  for  refined  meshes,  and  all  F,  g^,  /,  and  / .  an  estimate 

iiv'tmk  +  £  EM., .->>>] 

*  Arfll/IU  +  H/1k[o.rj  +  t  E  IMt, <».*<]  +  t  Ht*-') 

i*l  ^*8  t=l 

holds. 

Here  we  have  not  defined  all  our  norms,  and  have  altered  our  notation  [5] 
from  Section  2.  So  we  shall  explain  the  meaning  of  the  terms. 

The  first  term  on  the  left  is  the  discrete  iz(z)  norm  (just  defined)  of  the 
approximate  solution  at  time  t  =  T.  The  second  term  is  the  sum  of  discrete 
l2(t)  norms  of  the  approximate  solution  at  the  boundaries,  one  for  each  strip  St. 
The  first  term  on  the  right  is  the  discrete  f2(z)  norm  of  the  initial  function.  The 
second  term  is  the  discrete  l2(x,  t )  norm  of  the  inhomogeneous  forcing  function 
in  (2.11).  The  third  term  is  the  sum  of  discrete  l2(t)  norms  of  the  inhomogene¬ 
ous  boundary  terms  (2.12),  one  for  each  strip  5t.  Finally,  the  fourth  term  is  the 
sum  of  errors  which  arise  when  we  interpolate  from  a  coarse  mesh  to  a  finer  one 
during  mesh  adjustment  or  creation. 

For  a  uniform  mesh,  we  believe  that  a  scheme  is  GKS  stable  if  it  is  stable  in 
the  sense  of  Definition  3.6,  either  for  a  quarter-plane  or  strip  problem.  We 
believe  that  the  converse  is  not  true  in  general;  that  is,  the  new  definition  is 
stronger. 
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A  necessary  part  of  any  stability  proof  is  stability  of  two  quarter-plane  prob¬ 
lems  joined  along  a  coarse/fine  interface.  This  question  has  been  examined 
(using  the  GKS  definition)  in  Ciment  [9]  and  Oliger  [35]  for  the  case  of  equal  time 
steps  on  both  sides  of  the  interface.  Oliger  found  that  if  leap-frog  was  used  on 
both  sides  of  the  interface,  certain  restrictions  on  the  refinement  ratio  needed 
to  be  made.  But  if  the  difference  scheme  was  dissipative  on  one  side  of  the 
interface,  all  stability  problems  vanished.  Berger  [2]  proved  the  GKS  stability  of 
the  Lax-Wendroff  method  on  a  model  problem  along  a  coarse /fine  interface  with 
unequal  time  steps.  For  these  reasons  we  have  used  a  dissipative  scheme  in  all 
our  calculations. 

5.  Error  Analysis.  This  section  summarizes  partial  convergence  results 
which  give  a  theoretical  justification  for  our  algorithm. 

It  is  clear  from  the  description  of  our  algorithm  (and  from  the  computa¬ 
tions  in  Section  8.5)  that  mesh  refinement  does  not  increase  the  (global)  order 
of  accuracy  of  a  difference  approximation  (compared  to  using  a  similar  approxi¬ 
mation  on  a  uniform  mesh).  How  then  does  mesh  refinement  achieve  its 
efficiencies? 

In  [5]  we  gave  the  answer  using  a  proposition  which  we  did  not  prove,  but 
which  is  well-supported  by  computational  evidence  (see  Section  8).  It  is  the 
direct  analogue  of  Gustafsson’s  [21]  result  on  the  rate  of  convergence  of 
difference  approximations  to  the  initial  boundary  value  problem  for  hyperbolic 
systems  in  one  space  dimension.  Our  proposition  states  that  a  result  similar  to 
Gustafsson's  holds  when  (a)  we  replace  a  uniform  grid  by  our  refined  grid;  and 
(b)  the  difference  scheme  is  stable  according  to  Oliger’s  stability  definition 
rather  than  the  GKS  definition.  Using  the  norms  given  in  the  previous  section, 
this  proposition  states  that  the  discrete  lz  norm  of  the  global  error  at  time 
t  =  T  is  bounded  by  a  constant  Kf  times  a  sum  of  four  terms  These  terms  are 
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the  discrete  i2  norms  of  the  interpolation  error,  and  of  the  local  errors  of  the 
interior,  interface,  and  boundary  approximations. 

This  proposition  answers  the  question.  How  does  the  order  of  accuracy  of 
the  interpolation,  and  of  the  interior,  boundary  and  interface  approximations 
affect  the  global  error?  In  particular,  suppose  that  one  uses  an  0(hz )  approxi¬ 
mation  in  the  interior  of  refinements  and  at  coarse /fine  interfaces.  Also  sup¬ 
pose  that  one  uses  0(h)  boundary  approximations  and  linear  interpolation 
(0(/i2))  to  obtain  solution  values  on  f  +  l-st  level  refinements  from  f-th  level 
refinements.  Finally,  assume  that  this  scheme  is  stable  in  Oliger’s  sense.  Then, 
subject  to  certain  compatibility  and  smoothness  assumptions,  the  proposition 
states  that  the  global  error  is  0(hz). 

Now  suppose  we  use  a  strategy  suggested  by  de  Boor  [12]  and  Pereyra  and 
Sewell  [37]  in  other  contexts:  arrange  our  spatial  mesh  so  that  it  "approximately 
equidistributed"  the  hybrid  local  truncation  error  at  time  t  =  f<_l, 

i  =  1,2 . s;  compute  forward  in  time  until  the  equidistribution  condition  is 

"nearly  violated”  at  time  t  =  f*;  and  then  approximately  equidistribute  again. 
(The  hybrid  truncation  error  is  a  suitable  blending  of  the  interior,  interface  and 
boundary  truncation  errors.)  In  practice,  of  course,  it  is  more  work  to  discover 
whether  the  condition  is  "nearly  violated"  than  it  is  to  simply  approximately 
equidistribute  again.  That  is  why  we  choose  our  "equidistribution"  (mesh  adjust¬ 
ment)  times  a  priori.  A  somewhat  similar  strategy  has  been  used  by  Gannon 
[16]  for  parabolic  problems  with  finite  elements  in  two  space  dimensions. 

We  should  emphasize  that  we  do  not  approximately  equidistribute  in  prac¬ 
tice.  as  it  would  be  too  expensive,  Our  recursive  refinements  achieve  a  primitive 
form  of  approximate  equidistribution  at  much  less  cost. 

Using  our  proposition  and  results  of  Pereyra  and  Sewell  we  can  then  show 
(generalizing  what  Oliger  [33]  did  for  the  Cauchy  problem)  that,  loosely  speak- 
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ing,  using  our  mesh  refinement  algorithm  multiplies  the  constant  Kf  by  the  fac¬ 
tor 

Here 

K  = 

M=» 

is  an  upper  bound  on  the  ratio  maXj/i;/ minify  of  spatial  step  sizes  in  our 
refined  mesh;  is,  loosely  speaking,  the  ratio  of  the  length  of  the  spatial 
region  that  needs  refinement  (has  high  truncation  error)  to  the  length  of  the 
whole  spatial  region  b  -  a;  p  is  the  order  of  the  interior  approximation  and  the 
global  error;  and  q  and  r  are  the  stencil  sizes  given  in  (2.11).  When  the  solution 
has  rapid  variations  only  in  a  small  part  of  the  (spatial)  region,  then  the  local 
truncation  error  is  small  over  most  of  the  region,  and  fim „  is  therefore  small. 
Thus  our  algorithm  can  use  fewer  mesh  points  in  regions  where  the  local  trunca¬ 
tion  error  is  small  (compared  to  using  the  same  difference  scheme  on  a  uniform 
mesh  which  achieves  the  same  level  of  accuracy)  and  this  produces  significant 
economies,  as  shown  in  Section  8.4. 

6.  Data  Structures.  In  this  section  we  discuss  the  data  structures  used  in 
our  mesh  refinement  algorithm.  The  data  structure  has  two  parts:  a  four-way 
linked  tree  to  show  the  relationships  between  refinements,  and  sequentially  allo¬ 
cated  deques  for  storing  solution  values  of  refinements.  We  will  describe  the 
deques  first. 

6.1.  Deques.  One  of  the  most  important  operations  on  a  refinement  is  to 
"move"  it  left  or  right  to  follow  a  wave.  An  efficient  way  to  "move”  a  refinement 
right  (without  actually  moving  any  grid  points)  is  to  add  grid  points  to  the  right 
and  delete  them  from  the  left.  This  leads  us  to  store  the  solution  values  for  a 
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refinement  in  a  data  structure  called  a  "deque",  or  double-ended  queue  [28],  A 
deque  has  the  property  that  items  are  added  to  or  deleted  from  its  ends,  but  are 
never  inserted  in  or  deleted  from  the  middle. 

For  simplicity,  we  first  consider  a  scalar  equation.  A  natural  way  to  store  a 
collection  of  deques  is  sequentially,  as  shown  in  Fig.  4.  Here  we  see  a  region  with 
two  refinements,  and  one  of  the  refinements  itself  contains  a  refinement.  The 
solution  values  for  the  coarsest  mesh  occupy  a  fixed  region  at  the  lowest  end 'of 
a  vector  which  we  will  call  v .  The  solution  values  for  refinements  occupy  con¬ 
tiguous  sections  of  the  remaining  available  memory,  with  variable-width  gaps  of 
free  space  separating  them.  The  gaps  allow  us  to  expand  or  contract 
refinements  (to  a  limited  degree)  without  moving  the  solution  values.  The  solu¬ 
tion  values  corresponding  to  refinements  are  ordered  as  follows. 

The  coarse  mesh  is  labelled  refinement  1.  It  is  followed  in  v  by  the  "second 
level"  refinements  (labelled  2  and  3  in  Fig.  4),  which  arc  ordered  in  v  in  the 
same  order  as  the  refinements  are  encountered  in  proceeding  from  left  to  right 
in  the  computational  region.  Following  these  are  the  "third  level"  refinements 
(as  is  refinement  4  in  Fig.  4),  again  in  the  same  order  as  they  occur  in  a  left-to- 
right  scan  of  the  computational  region,  and  ignoring  the  positions  of  any  coarser 
(second  level)  refinements  encountered  in  the  computational  region.  Then 
would  appear  all  fourth  level  refinements,  and  so  forth.  This  scheme  duplicates 
certain  solution  values  in  the  vector  v,  namely  the  ones  which  correspond  to 
mesh  points  which  lie  on  different  level  refinements.  However,  doing  this  makes 
the  program  much  simpler. 

For  an  n  x  n  system  of  equations,  we  replace  the  solution  value  vector  v  by 
a  matrix  with  n  rows.  Each  row  has  the  same  arrangement  of  solution  values 
separated  by  gaps,  as  illustrated  in  Fig.  4  For  simplicity,  in  the  rest  of  this  sec¬ 
tion  we  shall  again  refer  to  v  as  a  vector. 
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6.2.  Hie  Tree.  Next  we  will  describe  the  (four-way  linked)  tree  of  records 
which  shows  the  relationships  between  refinements.  Trees  are  natural  here, 
since  we  use  recursive  refinements,  and  are  used  in  most  adaptive  solvers  for 
elliptic  equations.  There  is  a  one-to-one  correspondence  between  nodes 
(records)  of  the  tree  and  refinements,  with  the  root  corresponding  to  the  coar¬ 
sest  mesh.  In  the  following,  we  will  identify  a  refinement  with  its  node  (record), 
and  use  the  term  "refinement"  to  mean  "the  node  corresponding  to  a 
refinement".  We  will  sometimes  call  the  coarsest  mesh  a  "refinement"  for  uni¬ 
formity. 

The  root  of  the  tree  has  level  1,  its  immediate  successors  are  at  level  2, 
their  successors  have  level  3,  and  so  forth.  Each  node  contains  all  the  informa¬ 
tion  about  a  refinement,  except  its  solution  values.  We  will  now  describe  this 
information.  The  indices  base  and  top  indicate  where  in  the  vector  v  the  solu¬ 
tion  values  for  a  refinement  are  located.  This  is  shown  in  Fig.  4  for  the  fourth 
(level  3)  refinement,  but  omitted  for  the  other  refinements  to  avoid  clutter.  Also 
shown  is  a  pointer  coarse  to  the  parent  of  each  refinement.  Furthermore,  we 
need  pointers  to  all  refinements  ("children”)  of  a  refinement.  We  can  avoid  using 
a  variable  number  of  pointer  fields  for  this  by  using  the  usual  device.  We  use  one 
pointer  to  the  leftmost  descendant  (called  fine  in  Figure  6.1)  and  then  chain 
together  all  immediate  descendants  (children)  using  the  "right"  pointers,  called 
rlink  in  the  figure.  A  refinement  other  than  the  root  also  needs  two  indices  to 
denote  its  endpoints  within  its  parent,  that  is,  which  part  of  its  parent  it  refines. 
These  are  not  shown  in  the  figure. 

Since  we  will  often  add  or  delete  nodes  in  an  unpredictable  order,  we  imple¬ 
mented  the  tree  as  a  linked  list.  So  far  our  records  form  a  triply  linked  tree, 
exactly  as  in  Knuth  [28,  p.352].  However,  additional  links  are  needed. 


The  solution  is  advanced  in  time,  and  the  error  is  estimated  a  level  at  a 
time.  Because  we  already  have  the  rlink  pointers,  we  can  chain  together  all 
refinements  on  the  same  level  (not  just  those  with  a  common  parent)  using 
rlink.  Then  we  introduce  a  vector  of  pointers  pointing  to  the  leftmost 
refinement  on  each  level.  (These  are  shown  in  Fig.  4.)  This  is  related  to  the 
level-order  representation  of  a  tree  [28,  p.350]. 

The  last  operation  needed  on  our  data  structure  is  a  repacking  of  the  v  vec¬ 
tor.  to  be  discussed  shortly.  This  requires  us  to  sweep  through  v  in  both  direc¬ 
tions.  Thus  we  also  require  our  rlink  pointers  to  point  from  the  rightmost 
refinement  (node)  on  level  l  to  the  leftmost  refinement  on  level  l  +  1.  To  enable  a 
leftward  sweep,  we  introduce  "left"  pointers  Uink,  which  are  inverse  to  the  rlink 
pointers.  That  is,  if  node  p  has  right  pointer  rlink  pointing  to  node  q ,  then  q  has 
left  pointer  llink  pointing  to  p . 

The  final  result  is  a  four-way  linked  tree:  a  triply  linked  tree  with  the  addi¬ 
tional  property  that  all  the  nodes  are  linked  together,  in  level  order,  in  a  doubly 
linked  list.  Our  tree  uses  an  inconsequential  amount  of  memory,  compared  to 
the  memory  ini'  devoted  to  solution  values. 

We  now  examine  how  the  operations  on  refinements  are  effected  using  this 
data  structure.  Advancing  the  difference  approximation  (in  time)  or  estimating 
the  error  can  be  done  a  level  at  a  time,  using  the  rlink  pointers  and  the  leftmost 
pointers  on  each  level.  Here  we  also  use  the  "ancestor"  or  coarse  pointers  to 
copy  solution  values  from  finer  meshes  to  coarser  meshes  for  points  *  which  lie 
on  more  than  one  refinement. 

Similarly,  we  adjust  the  refinements  level  by  level,  starting  with  the  highest 
(finest)  level.  The  mesh  adjustment  operations  can  be  effected  using  four  ele¬ 
mentary  operations,  which  are  natural  for  a  deque.  They  are  shorten  left,  shor¬ 
ten  right,  extend  left,  and  extend  right.  Shortening  either  end  of  a  refinement  is 
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a  trivial  operation,  accomplished  by  moving  a  base  or  top  index.  Deleting  a 
refinement  is  the  sarnie,  but  also  involves  removing  a  record  from  the  tree  and 
reclaiming  its  storage.  If  there  is  enough  space  available,  extending  either  end 
of  a  refinement  involves  changing  an  index,  copying  solution  values  from  the 
parent  refinement,  and  filling  in  new  solution  values  using  linear  or  quadratic 
interpolation  in  space.  Creation  is  the  same,  plus  the  operation  of  inserting  a 
new  node  in  the  tree.  Separation  of  a  refinement  into  two  refinements  involves 
changing  indices  and  inserting  a  new  node.  Finally,  merging  two  refinements  is 
simplified  because  we  insisted  on  the  left-to-right  ordering  of  refinements  in  the 
vector  v .  We  move  left  the  solution  values  of  the  right  refinement,  if  necessary, 
then  extend  the  left  refinement  to  the  right,  change  some  indices,  and  delete 
the  right  node.  Complicating  the  last  two  operations  is  the  need  to  adjust 
pointers  to  descendant  refinements. 

6.3.  Memory  Repacking.  A  problem  occurs  during  an  "extend”  operation 
when  there  is  insufficient  expansion  room  between  refinements.  This  cetlls  for  a 
repacking  of  memory,  and  an  edgorithm  for  this  is  given  by  Knuth  [28,  pp. 245-6] 
for  the  case  of  a  sequence  of  stacks  (rather  than  deques).  We  will  therefore 
describe  the  modifications  to  this  algorithm  for  our  data  structure. 

When  a  refinement  runs  out  of  room  in  the  v  vector,  moving  only  the  adja¬ 
cent  refinement  will  probably  cause  another  repacking  to  occur  soon,  so  it  is 
better  to  reallocate  all  available  memory  when  a  refinement  runs  out  of  room. 
Knuth  breaks  this  into  two  parts:  Algorithm  G,  which  decides  how  to  allocate  the 
free  memory  to  the  refinements,  and  Algorithm  R,  which  actually  moves  the 
refinements  into  the  positions  dictated  by  Algorithm  G.  It  is  Algorithm  R  which 
requires  the  forward  and  backward  sweep  of  the  v  vector  in  order  to  avoid 
overwriting  any  information. 
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We  used  Algorithm  R  unchanged  and  modified  Algorithm  G  as  follows 
Khuth's  maun  idea  is  to  share  ten  percent  of  the  free  memory  equally  among  the 
refinements,  and  to  divide  the  other  ninety  percent  proportionate  to  the  amount 
of  increase  in  refinement  size  since  the  previous  repacking.  This  idea  is  not  use¬ 
ful  in  our  case.  For  a  traveling  wave,  all  refinements  stay  about  the  same  size, 
but  "move".  However,  we  can  modify  this  rule  by  awarding  the  ninety  percent  of 
available  memory  proportionate  to  the  amount  each  refinement  has  moved 
since  the  last  repacking.  We  discover  whether  a  refinement  has  moved  primarily 
left  or  right  (in  memory)  since  the  last  repacking,  and  award  its  share  of  the 
ninety  percent  to  its  left  or  right,  respectively. 

This  change  to  Knuth's  algorithm  greatly  reduces  the  number  of  repackings 
compared  to  more  naive  allocation  methods.  Since  the  coarse  mesh  doesn’t 
move  it  receives  none  of  the  ninety  percent  allocation  Furthermore,  the  higher 
level  refinements  move  further  (measured  in  number  of  mesh  points,  not  physi¬ 
cal  distance)  than  the  lower  level  ones,  so  they  are  awarded  more  free  space  by 
this  scheme. 

Our  storage  scheme  for  solution  values  cannot  be  generalized  to  more 
space  dimensions,  because  refinements  no  longer  have  only  two  "ends".  In  more 
dimensions  it  is  also  unwise  to  merge  adjacent  refinements,  even  if  they  are  on 
the  same  level.  But  our  scheme  avoids  more  complicated  storage  allocation 
algorithms  [2]. 

7.  Choice  of  Programming  Language.  We  will  now  explain  and  justify  our 
choice  of  implementation  method  for  the  programs  used  in  our  computations. 

Possible  alternatives  include  Algol  W.  PL/I,  Algol  60,  Algol  68.  Pascal.  For¬ 
tran,  and  Fortran  with  preprocessor.  The  arguments  against  the  first  four  are 
lack  of  availability  of  a  compiler  and/or  lack  of  portability  Raw  Fortran  (even 
Fortran  77)  is  cumbersome  to  use  because  of  the  lack  of  control  and  data 
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structures,  both  of  which  are  crucial  for  our  task.  However,  most  numerical 
software  is  written  in  Fortran,  and  if  one  uses  another  language,  there  must  be 
an  interface  to  Fortran.  For  reasons  given  in  Kernighan  [27],  we  did  not  use  Pas¬ 
cal,  despite  its  excellent  data  structuring  facilities. 

We  then  examined  preprocessors.  We  rejected  Feldman’s  [15]  EFL,  Grosse’s 
[20]  language  T,  and  Stein's  [45]  language  for  portability  reasons,  even  though 
their  respective  authors  have  devoted  considerable  thought  to  devising 
appropriate  language  constructs  for  numerical  algorithms.  We  examined  two 
portable  Fortran  preprocessors:  Kernighan's  [26]  Ratfor  and  Cook  and  Shustek’s 
[10]  Mortran.  Although  the  former  is  more  widely  used,  we  chose  the  latter 
because  it  is  far  more  general  and  flexible.  (Engquist  and  Smedsaas  [14]  have 
also  used  a  macro  processor  in  their  work.  Gropp  [19]  has  developed  a  language 
for  mesh  refinement  algorithms.) 

The  term  Mortran.  like  Fortran,  has  several  meanings.  It  can  mean  a  struc¬ 
tured  source  language,  a  translator  for  that  language,  or  a  macro-processor. 
The  structured  language  is  implemented  as  a  set  of  macros  which  are  used  by 
the  Mortran  macro  processor  to  translate  the  language  into  Fortran.  The  result¬ 
ing  Fortran  program  is  then  run  like  any  other  Fortran  program. 

In  contrast  to  most  other  Fortran  preprocessors,  the  Mortran  preprocessor 
is  written  in  a  portable  subset  of  ANSI  (standard)  Fortran.  Hence  the  Mortran 
preprocessor,  and,  more  importantly,  Mortran  source  programs,  are  portable 
between  different  machines.  (We  ran  our  programs  on  an  IBM  370/168,  a  CDC 
7600,  and  a  DEC  VAX  with  minimal  conversion  problems.)  Furthermore,  Mortran 
source  and  Fortran  source  can  be  intermixed,  so  the  Mortran  user  has  access  to 
all  existing  Fortran  software. 

We  felt  that  there  was  one  property  of  Mortran  which  made  it  especially 
desirable  for  this  project:  extensibility.  This  means  that  new  data  structures, 
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operations  on  data  structures,  and  control  structures  can  be  added  to  the 
language  (at  rather  small  cost  in  implementation  time)  by  adding  additional 
macros  to  the  language.  To  implement  the  linked  list  for  our  four-way  linked 
tree,  we  needed  records  and  pointers,  and  Mortran  allows  us  to  create  these  new 
data  types  and  to  define  operations  (such  as  following  pointers)  on  these  data 
types.  We  modified  the  Mortran  macros  for  records  and  pointers  given  in  Zahn 
[49],  and  used  Pascal-like  syntax  [25]. 

Of  course,  many  of  the  restrictions  of  Fortran  remain.  Among  these  are 
lack  of  dynamic  storage  allocation,  arrays  with  arbitrary  subscripts,  and  recur¬ 
sion.  Except  for  recursion,  these  are  not  serious.  We  needed  recursion  when 
advancing  the  solution  and  when  plotting  it;  for  the  latter  we  needed  to  search 
the  tree  in  preorder  to  produce  solution  values  ordered  with  increasing  *.  Doing 
without  recursion  made  for  more  obscure  code. 

For  most  control  constructs.  Mortran  produces  Fortran  code  that  is  as 
efficient  as  possible  without  using  global  flow  analysis.  For  one  type  of  loop  (the 
for  loop)  we  rewrote  the  macros  to  generate  more  efficient  Fortran  [5]. 

A  criticism  often  levelled  at  macro  preprocessors  is  that  they  produce 
error  diagnostics  in  terms  of  Fortran  instead  of  the  source  language.  Zahn  [49] 
shows  how  to  write  additional  macros  to  produce  reasonable  source-level  diag¬ 
nostics  for  the  record  and  pointer  constructs  These  macros  also  insure  that  a 
pointer  points  only  to  a  record  of  the  appropriate  type. 

In  [5]  we  gave  a  complete  program  listing  of  our  mesh  refinement  algorithm 
applied  to  problem  P2  in  Section  6. 

6.  Computational  Results.  In  this  section  we  answer  the  following  questions 
about  our  method: 
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1.  Does  our  method  "follow"  or  "track"  steep  gradients?  Is  it  fooled  by  back¬ 
ground  effects? 

2.  Is  the  algorithm  sufficiently  general  to  allow  refinements  to  be  created,  des¬ 
troyed,  merged,  separated,  moved,  and  to  abut  boundaries? 

3.  Is  the  algorithm  sensitive  to  the  direction  of  characteristics,  or  dependent 
on  knowing  that  certain  boundary  conditions  are  inflow  or  outflow? 

4.  Will  the  method  handle  nonlinear  problems? 

5.  How  well  will  the  method  follow  discontinuities  or  shocks? 

6.  Are  recursive  refinements  worthwhile? 

7.  How  should  one  choose  the  refinement  ratios  N  and  Ml 

8.  How  efficient  is  the  method,  both  in  execution  time  and  memory? 

9.  How  does  the  global  error  behave  as  h  -»  0? 

10.  How  do  the  two  methods  of  interior  local  error  estimation  of  Section  3  com¬ 
pare  in  accuracy  and  efficiency? 

11  How  do  different  boundary  approximations  and  methods  of  estimating  their 
error  affect  the  solution? 

12.  How  often  should  one  monitor  the  local  truncation  error  (and  adjust 
refinements)? 

8.1  Model  Problems.  Since  we  believe  it  is  impossible  to  answer  these  ques¬ 
tions  analytically,  we  resort  to  numerical  experiments  on  model  problems. 
Problem  Pi,  the  first-order  wave  equation,  was  introduced  in  Section  2.4.  We 
now  introduce  two  additional  problems. 

P2  is  the  second  order  wave  equation,  written  as  a  2  x  2  first  order  system, 
with  "open”  boundary  conditions  (i.e. ,  the  boundaries  are  "transparent"  to  trav¬ 
eling  waves)  As  exact  solution  wc  use  two  counter-streaming  Gaussian  pulses, 


superimposed  on  a  sinusoidal  background.  The  differential  equation  is 


ut  =  AUg,  a  s  x  s  6 , 0  £  t ,  0  <  c  , 

where 

fo  c 

*  =  [c  OJ* 

with  initial  conditions 


ui(*.0)  =  f(x)+g(x) 
u2(x,0 )  =  -/(x)  +  $(x) 


and  open  boundary  conditions 


U](a,f)  =  u2(a,f)  +  2/(a  -  cf) 
u^b.t)  -  —u2(6,f)  +  2 g{b  +  ct) 


12  0. 


(PI) 


We  choose  a  =  0,  6  *  4,  c  =  1.  The  exact  solution  is 

Uj(x,f)  =  /(x  -  cf)  +  g{x  +  cf), 
u2(x,f)  =  -/(x  -  ct )  +  £ (x  +  cf ). 

To  produce  our  interacting  pulses,  we  take  /  (x)  exactly  as  in  PI.  and 
g{x)  =  -exp(~a(x-4.5)2), 

where  a  =  200.  Each  pulse  occupies  about  8  percent  of  the  region  a  «  x  ss  6 . 

The  difference  approximation  in  the  interior  is  again  Lax-Wendroff 
i^(f+Jbt)  =  (/  +  AktD\>  +  KAWD^Dityit) 

(omitting  superscripts  f  onv),  with  coarse /fine  approximation 
(*+*,)  =  (/■*■  AkLD\>-'  + 

at  interfaces,  the  obvious  initial  condition,  and  obvious  boundary  approxima¬ 
tions  for  V).  For  v2,  we  need  extra  boundary  conditions  at  both  x  =  o  and  x  =  6 , 


49 


and  we  use  either 

(a)  upwind /downwind  differencing: 

Vj2(t+kt)  x  (/  +  cktD+)Vjz(t)  at  x  x  a, 
v}-2(t+ki )  =  (I  +  ckt Di )vi2(t )  at  x  =  6. 

where  viz (<)  denotes  an  approximation  to  uz(x  ,t)  at  x  =  a  +  fa  on  an  i-th  level 
refinement;  or 

(b)  first-order  extrapolation 

(D+fojzit  +kt)  =  0  at  x  x  a. 

{Dl-Yvjz{t  +kt)  =  0  at  x  x  b. 

(The  first  2  in  each  line  is  an  exponent,  not  a  superscript.)  Gustafsson.  Kreiss 
and  Sundstrom  [22]  showed  that  both  approximations  (a)  and  (b)  are  stable  with 
Lax-Wendroff,  according  to  their  stability  definition. 

Our  second  problem  is  the  Riemann  shock  tube  problem  for  a  compressible, 
inviscid  gas  [43].  The  compressible  Euler  equations,  in  conservation  form,  are 

Ut  +  F(U)X  =0,  OszsUaO, 

where 


p 

m, 

u  = 

m 

.  and  F(U)  = 

(mz/p)  +  p 

e 

(e  +p)m/p 

Here  p.  m,  and  e  are  the  mass,  momentum,  and  energy,  respectively,  per  unit 
length,  andp  is  the  pressure.  Let  u  =  m/ p  be  the  velocity.  If  the  gas  is  polytro¬ 
pic  with  ratio  of  specific  heats  y,  we  can  express  the  energy  as  e  =  p(e  +  )fcuz), 
where  the  internal  energy  per  unit  mass  t  -p/(y-  1  )p.  We  can  also  write 
P  =  [y  ~  1)(«  -)$mz/p). 
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Initially  we  will  have  two  gases  at  rest,  but  at  different  pressures  and  densi¬ 
ties,  separated  by  a  membrane.  At  time  t  =  0  the  membrane  is  instantaneously 
removed.  Following  Sod  [43],  let  us  define  the  left  state  S{  =  (1.0,  0..  2.5)r  and 
the  rigkJ.  state  Sr  =  (.125,  0.,  .25)r.  Then,  initially,  U(x,  0)  =  St  for  0  x  <  0.5, 
and  U(x,  0)  =  5r  for  0,5sxs  1.  At  the  left  boundary  U(0,  t)  =  St,  and  at  the 
right  boundary  U(l,  t)  -  Sr,  both  for  t  >  0.  but  before  reflections  off  the  wall. 
(One  should  actually  specify  only  one  of  p,  m,  or  e  at  each  boundary,  as  can  be 
seen  by  considering  characteristics  and  Riemann  invariants.  However,  our  pro¬ 
cedure  is  equivalent  before  reflections  off  boundaries,  since  the  solutions  are 
constant  near  the  walls.)  For  fixed  t  >  0,  before  reflections  off  the  walls,  the 
exact  solution  is  a  function  of  (x  -  0.5)/ 1  only.  Proceeding  from  left  to  right 
across  the  region,  the  exact  solution  is  the  constant  state  St,  followed  by  a  rare¬ 
faction,  followed  by  a  constant  state.  To  the  right  of  this  is  a  contact  discon¬ 
tinuity.  followed  by  another  constant  state,  then  a  shock,  concluded  by  the  con¬ 
stant  state  Sr.  Sod  shows  how  to  compute  the  exact  solution  at  any  time  by 
solving  a  3  by  3  system  of  nonlinear  equations. 

We  approximated  this  using  the  usual  two-step  Lax-Wendroff  method,  as 
given  in  Richtmyer  and  Morton  [41,  p.302],  We  added  dissipation  by  subtracting 
from  the  flux  occurring  in  the  two-step  Lax  Wendroff  scheme  the  quantity 
yfuj*+i  -  u/1(UJ*+,  -  U”),  with  v  -  1.  That  is,  the  dissipation  is  based  entirely  on 
the  values  of  U  and  u  at  the  old  time  level.  We  chose  this  method  rather  than 
Sod’s  because  it  is  easier  to  implement.  At  the  initial  jump  discontinuity  we 
took  average  values  for  the  density  and  energy.  Naturally  we  modified  all  this  by 
adding  subscripts  and  superscripts  l  (for  level  of  refinement)  in  the  appropriate 
places.  Since  our  previous  calculations  used  time-dependent  boundary  condi¬ 
tions,  we  use  constant  ones  here.  The  solution  was  calculated  up  to  and  includ¬ 
ing  time  t  =  0.15.  (Of  course,  at  coarse/fine  interfaces,  we  made  modifications 
as  in  Pi  and  P2.) 
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8.2  Qualitative  Results.  Figs.  5(a)  through  5(o)  show  our  algorithm  applied 
to  the  second-order  wave  equation  with  counter-streaming  pulses.  We  used  a 
coarse  mesh  of  81  points,  with  refinement  ratios  N  =  M  -  3  and  maximum 
number  of  refinement  levels  =  5.  (Only  four  levels  were  actually  used  in  the  cal¬ 
culation.)  The  local  error  tolerance  was  0.01,  Aj  =  0.8,  and  the  wave  speed  c  -  1. 
We  used  the  Richardson  method  to  estimate  the  error. 

Each  graph  plots  the  approximate  solution  component  t \(z,t),  i  -  1,  2, 
versus  x  (at  a  fixed  time)  as  a  solid  line,  together  with  a  separate  calculation 
done  with  no  refinement  (maximum  refinement  level  =  1),  shown  as  a  dotted 
line.  We  usually  show  the  first  component  v^z.t)  versus  x,  but  in  a  few 
instances  we  show  vz(x,t)  versus  x.  Since  the  maximum  error  in  this  calcula¬ 
tion  was  less  than  2  percent,  the  refined  solution  may  be  taken  as  the  exact 
solution  on  the  graph. 

In  Fig.  5(a)  the  pulses  have  not  yet  entered  the  region  and  we  see  only  the 
sinusoidal  background.  (This  graph  was  made  at  f  =  0.04  rather  than  f  =  0  since 
we  used  the  exact  solution  at  t  =  0.04  to  compare  with  another  method  not 
given  here.  We  started  at  t  =  0  in  problem  P3.)  In  Fig.  5(b)  both  pulses  have 
entered  and  refinements  on  levels  2,  3  and  4  have  been  created  at  both  boun¬ 
daries.  (The  small  numbers  at  the  top  are  the  level  numbers  of  the  ends  of 
refinements.)  All  the  refinements  abut  the  boundaries.  In  Fig.  5(c)  both  pulses 
have  left  the  boundaries  and  are  moving  towards  each  other.  The  refinements 
follow. 

In  Fig.  5(d)  the  second-level  refinements  are  about  to  merge,  and  in  Fig. 
5(e)  they  have  merged;  however,  the  third  and  fourth  level  refinements  have  not 
yet  merged  Note  that  the  unrefined  solution ' ’s  become  a  very  poor  approxi¬ 
mation  to  the  pulses— the  unrefined  peak  and  trough  have  only  half  the  size  they 
should.  Note  also  that  behind  the  pulses,  the  unrefined  solution  has  a  large 
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undershoot.  In  Fig.  5(f)  the  third  level  refinements  have  merged.  Fig.  5(g) 
shows  the  second  component  of  the  solution  at  the  same  time.  In  Fig.  5(h)  the 
fourth  level  refinements  have  merged.  Fig.  5(i)  shows  the  second  component  of 
the  solution  at  the  same  time. 

In  Fig.  5(j)  the  pulses  have  crossed,  and  the  third  and  fourth  level 
refinements  have  separated.  (Note  that  the  pulses  cross,  but  the  refinements  do 
not.)  In  Fig.  5(k)  the  second  level  refinements  have  separated  as  well.  Note  the 
degradation  in  the  unrefined  solution  at  this  point.  Fig.  5(1)  shows  the  pulses 
approaching  the  boundaries.  Now  the  unrefined  solution  has  phase  errors  as 
well  as  amplitude  errors.  In  Fig.  5(m)  the  pulses  are  leaving  the  region,  and  he 
refinements  again  abut  the  boundary.  Fig.  5(n)  shows  that  the  fourth  level 
refinements  have  been  deleted.  Finally,  Fig.  5(o)  shows  that  all  refinements  have 
been  deleted.  Only  the  sinusoidal  background  remains. 

This  problem  answers  questions  1  to  3.  It  show's  the  interactions  of  up  to 

seven  refinements.  Since  both  boundaries  act  as  inflow  boundaries  at  some 

times  and  as  outflow  boundaries  later,  our  method  does  not  depend  on  the  direc¬ 
tion  of  characteristics.  These  calculations  also  show  that  the  method  is  not 
fooled  by  background  oscillations,  and  that  the  refinements  "follow”  and  resolve 
steep  gradients.  The  method  clearly  adapts  to  time-dependent  boundary  condi¬ 
tions. 

Figs.  6(a)  to  6(d)  show  the  algorithm  applied  to  the  Riemann  shock  tube 
problem  (P3).  Shown  at  time  t  =0.15  are  the  density,  velocity,  pressure  and 
internal  energy  vs.  x.  The  dotted  line  shows  the  unrefined  calculation.  We  used 
101  points  on  the  coarsest  mesh,  with  refinement  ratios  N  =  M  -  4,  and  max¬ 
imum  number  of  levels  =  3.  The  local  error  tolerance  was  0.01,  and  the  max¬ 
imum  wave  speed  input  to  the  program  was  max  ((u|  +  c)  a  2.2.  We  took 

X,  =  15/37  «  .405405,  since  X  =  5/  11  »  .45454  is  necessary  for  CFL  stability  in 


the  absence  of  dissipation.  The  maximum  number  of  reflnment  levels  was  3.  We 
used  the  Richardson  method.  (In  contrast  to  the  smooth  solution  case,  the  local 
truncation  error  estimate  around  a  discontinuity  does  not  decrease  on  finer 
meshes;  hence  our  method  will  always  use  the  maximum  number  of  levels  in  this 
case.)  The  purpose  of  this  calculation  is  to  see  if  our  method  can  handle  a  non¬ 
linear  problem,  and  if  it  can  '‘follow"  and  resolve  shocks,  contact  discontinuities, 
and  rarefactions,  for  which  it  was  not  designed.  These  results  are  directly  com¬ 
parable  with  those  of  Sod  [43]. 

The  shock,  rarefaction  and  sonic  point  (at  the  right  side  of  the  rarefaction) 
are  very  well  resolved.  The  contact  discontinuity  has  wiggles,  since  no  dissipa¬ 
tion  was  added  around  it.  The  maximum  error  occurs  here  also.  In  this  problem 
most  of  the  error  comes  from  the  initial  step  function,  and  our  method 
"instantly"  refines  around  the  initial  step.  The  method  has  placed  the  third  level 
refinements  only  where  needed  -  around  the  shock,  contact,  and  the  edges  of 
the  rarefaction. 

This  problem  shows  that  the  Richardson  error  estimation  method  performs 
correctly  under  far  more  general  circumstances  than  the  hypotheses  of 
Theorem  1.  We  need  not  assume  constant  coefficients,  linearity,  or  a  Cauchy 
problem.  This  method  performs  well  even  when  the  global  error  has  a  lower 
order  than  the  local  error,  and  when  the  solution  is  not  continuous.  (It  is  well- 
known  [31]  that  the  global  error  is  not  second  order  in  those  parts  of  the  region 
reachable  by  characteristics  emanating  from  the  shock.)  In  fact,  for  smooth 
solutions  this  method  yields  estimates  with  approximately  10  to  15  percent 
error.  (One  of  the  methods  examined  in  [5]  gave  800  to  1000  percent  errors!) 

Of  course,  we  cannot  yet  seriously  suggest  this  as  a  method  for  computing 

shocks,  since  we  have  not  considered  matters  such  as  entropy,  monotonicity,  or 

\ 

conservation.  (But  see  [2]  for  a  discussion  of  conservation  )  This  problem  is  too 
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simple  a  model  for  problems  with  discontinuous  solutions,  since  there  are  no 
collisions  of  shocks  and/or  rarefactions.  For  some  problems  (e.g.,  chemical 
kinetics)  it  is  important  to  accurately  resolve  shocks.  For  others  (e.g.,  gas 
dynamics)  our  method  (using  its  present  error  estimation)  may  expend  too 
much  effort  resolving  shocks.  In  the  latter  case  a  method  such  as  that  of  Wood¬ 
ward  and  Colella  [48]  may  be  more  appropriate.  One  could  also  construct  a 
hybrid  scheme,  by  applying  a  difference  method  designed  for  shocks  only  on 
high  level  refinements  containing  the  shock(s),  and  using  a  more  inexpensive 
difference  method  on  the  rest  of  the  refinements. 

We  now  proceed  to  more  quantitative  questions. 

8.3  Choosing  Refinement  Ratios  and  Maximum  Level.  In  [5]  we  studied  the 
questions  of  how  to  choose  the  refinement  ratios  N  and  M  and  whether  to  use 
recursive  refinements  with  computations  on  model  problem  PI.  We  found  that: 

(a)  It  is  necessary  to  use  recursive  refinements  for  efficiency. 

(b)  If  we  fix  the  coarse  mesh  size  and  the  local  error  tolerance,  and  let  the 
maximum  number  of  refinement  levels  increase,  then,  for  smooth  solutions,  the 
number  of  refinement  levels  used  first  increases  and  then  stays  constant.  That 
is.  the  method  only  uses  as  many  levels  as  are  necessary.  For  smooth  solutions, 
one  should  therefore  set  the  maximum  refinement  level  at  some  large  number. 
(In  many  cases,  the  method  uses  three  levels.) 

(c)  One  can  choose  the  refinement  ratios  N  and  M  to  be  between  4  and  6  for 
near  optimal  efficiency.  Nothing  is  gained  by  choosing  different  refinement 
ratios  for  different  l,  or  by  choosing  N*M. 

We  next  come  to  the  most  important  result  of  this  paper. 

8.4  Efficiency  of  the  Method.  In  Section  8.2  we  showed  that  our  method  was 
able  to  resolve  steep  gr  dients,  and  even  shocks,  in  the  solution  We  showed  this 
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by  comparing  the  solution  obtained  by  our  method  with  one  obtained  on  a  uni¬ 
form  coarse  mesh.  This  clearly  showed  the  qualitative  superiority  of  our 
"refined"  solution  over  the  "unrefined"  one. 

However,  the  "unrefined"  solution  cost  far  less  to  compute  than  the 
"refined"  one.  For  example,  it  cost  1.17  seconds  to  compute  a  "refined"  solution 
of  Problem  Pi  up  to  time  t  =  3.6  (without  graphic  output,  etc.)  vs.  0.04  seconds 
to  compute  the  "unrefined”  solution  up  to  the  same  time.  This  is  a  factor  of  29 
more  expensive.  But  the  unrefined  solution  was  worthless. 

Thus,  to  study  the  efficiency  of  our  method,  we  need  to  compare  the  com¬ 
puting  time  taken  by  our  method  with  the  computing  time  taken  to  produce 
(approximately)  the  same  error  on  a  uniform  (fine)  mesh.  As  a  by-product,  we 
will  also  be  able  to  compare  the  memory  taken  in  the  two  approaches. 

Our  method  is  simple.  Instead  of  comparing  a  "refined”  solution  E  with  a 
solution  computed  on  only  its  coarsest  mesh,  we  compare  E  with  a  solution  com¬ 
puted  on  a  uniform  fine  mesh  whose  spacing  is  the  same  as  the  spacing  of  the 
E's  finest  mesh.  If  these  two  produce  approximately  the  same  error,  then  we 
have  a  valid  comparison. 

Because  this  is  probably  the  most  important  result  in  this  paper,  we  made 
this  study  on  all  three  of  our  model  problems.  The  result  is  approximately  the 
same  for  all. 

Table  1  shows  the  results.  The  errors  are  at  t  -  3.6  for  PI  and  P2,  and  at 
t  =0.15  for  P3.  For  Pi  and  P2,  this  is  the  time  just  before  the  pulse  leaves  the 
spatial  interval.  (We  will  measure  the  error  at  these  times  in  the  rest  of  this  sec¬ 
tion.)  The  local  error  tolerance  was  0.001  for  PI  and  P2  and  0.01  for  P3.  We  used 
Richardson  extrapolation.  The  maximum  error  is  the  maximum  over  all  com¬ 
ponents  of  the  (global)  error  n-vector,  and  over  all  grid  points  existing  at  a 
given  time.  The  Z2(x)  norm  of  v  was  given  in  Section  4.  Here  we  alter  the 


definition  by  taking  the  ls  norm  of  each  component  of  the  solution,  and  then  tak¬ 
ing  the  maximum  of  the  n  results.  The  i2(z)  norm  of  the  error  is  analogous.  We 
used  upwind /downwind  boundary  conditions  for  PI  and  P2  (see  Section  8.1). 

In  our  tables,  a  maximum  mesh  level  of  1  signifies  that  only  the  coarse 
mesh  is  present  (no  refinement),  2  signifies  one  additional  level  of  refinement, 
and  so  forth.  The  times  reported  are  CPU  seconds  on  a  CDC  7800.  (Since  this 
machine  runs  only  in  batch,  these  times  are  highly  reproducible  between  runs, 
and  we  did  not  need  to  take  average  times.)  The  storage  used  is  for  solution 
values  only,  and  is  the  maximum  storage  used  at  the  (refinement)  level  listed, 
per  solution  component,  for  ail  times.  Since  the  coarse  mesh  is  static,  it  always 
uses  61  or  101  locations.  The  total  listed  is  the  sum  over  all  levels. 

We  see  that  in  terms  of  computer  time  our  method  is  3  to  5.5  times  as 
efficient  as  using  a  uniform  fine  mesh  which  produces  the  same  error.  In  terms 
of  memory,  a  factor  of  1.7  to  2.2  is  gained.  At  first  it  might  seem  surprising  that 
our  method  could  be  more  efficient,  since  it  requires  much  greater  overhead  (to 
estimate  the  error  and  adjust  refinements)  than  the  uniform  mesh  method.  This 
is  compensated  for,  however,  by  being  able  to  take  large  time  steps  in  unrefined 
regions.  Most  of  the  computing  time  is  spent  in  advancing  the  solution  on  the 
finest  refinements.  The  error  estimation  and  mesh  manipulations  are  almost 
free.  That  is  why  it  is  so  important  to  refine  only  where  necessary. 

In  general,  of  course,  the  specific  efficiency  factor  depends  primarily  on  the 
fraction  of  the  region  needing  refinement,  and  other  factors  such  as  the  local 
error  tolerance,  when  (for  which  t)  we  are  doing  the  comparison,  the  wave 
speed,  and  so  forth. 

These  results  show  why  we  need  to  create,  delete,  merge,  and  separate 
refinements.  If  we  did  not  allow  merging  and  separating,  then  in  P2,  we  would 
have  had  to  refine  the  whole  region  when  the  pulses  were  entering  or  leaving  the 
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region  and  this  would  degrade  the  efficiency.  Similarly,  in  P3.  we  would  have  had 
to  use  a  single  refinement  encompassing  all  the  non-constant  features  of  the 
solution. 

Our  mesh  refinement  algorithm  reduces  the  (maximum)  number  of  mesh 
points  needed  to  achieve  a  given  accuracy,  and  this  naturally  reduces  the 
amount  of  work,  but  does  the  amount  of  effort  per  mesh  point  decrease?  Table  1 
also  gives  this  figure,  obtained  by  dividing  the  computer  time  by  the  maximum 
number  of  mesh  points  used.  It  is'  clear  that  in  all  cases  the  work  per  mesh 
point  is  decreased  by  a  factor  of  approximately  two  (the  notation  n-m  means 
n*10m). 

It  might  be  argued  that  we  obtained  our  results  only  by  adjusting  or  tuning 
the  parameters  N,  M.  To  refute  this  charge,  we  have  shown  several  different 
values  of  N  and  M.  This  shows  that  although  we  cannot  easily  determine  the 
optimal  N,  M,  even  suboptimal  choices  still  yield  a  significant  savings  in  execu¬ 
tion  time. 

8.5  Behavior  as  h  -*  0.  For  our  mesh  refinement  algorithm  we  can  study 
two  types  of  convergence.  In  all  cases  we  let  A  =  A:t/At  =  constant,  independent 
of  t .  Thus  N  and  M  are  fixed.  We  shall  assume  the  exact  solution  is  sufficiently 
smooth. 

(a)  We  can  keep  the  local  error  tolerance  6  and  the  maximum  refinement 
level  A  fixed,  and  let  the  largest  spatial  step  h  *  A,  approach  zero.  If  we  take  a 
sufficiently  large  value  for  A.  then  the  algorithm  will  refine  as  much  as  it  needs 
to.  Furthermore,  (for  smooth  solutions)  our  method  then  has  a  property  which 
leads  to  simplified  analysis1.  For  sujyiciently  small  A,,  our  refined,  mesh  becomes 
a  uniform  mesh.  This  type  of  convergence  is  not  desirable,  since  the  advan¬ 
tages  of  refinement  are  ultimately  lost.  ' 


(b)  In  the  second  method  we  let  h-i  -*  0  and  choose  6  as  a  function  of  A,,  so 
that  6  -*  0  also.  (Alternatively,  we  could  let  6  -*  0  and  choose  A  j  as  a  function  of 
<5.)  If  one  knows  the  order  of  the  global  error,  one  could  choose  6  =  C(A1)p  for 
some  constant  C.  This  is  certainly  the  theoretically  most  appealing  method.  If 
we  use  this  method,  then  the  grid  does  not  approach  a  uniform  coarse  grid  as 
At  -»  0.  Rather,  the  ratio  of  the  width  of  any  refinement  to  the  width  of  its 
parent  should  approach  a  constant  as  Aj  -*  0.  However,  for  checking  the  asymp¬ 
totic  behavior  of  the  program  we  shall  first  use  method  (a),  since  it  does  not  beg 
the  question  by  assuming  the  behavior  of  the  global  error. 

We  first  keep  the  maximum  number  of  refinement  levels  constant,  and  less 
them  necessary  (for  the  method  to  refine  as  much  as  possible),  and  study  the 
first  type  of  convergence.  Table  2  shows  these  results  on  Pi  using  =  0.8, 
./V  =  M  =  4,  three-level  Richardson  extrapolation,  and  local  error  tolerance  = 
0.001.  For  the  smaller  values  of  A,  the  0(h2)  behavior  of  the  errors  (both  max¬ 
imum  and  l2)  is  apparent. 

Next  we  do  the  same  test,  but  choose  the  maximum  number  of  levels  large 
enough  so  that  the  method  refines  as  much  as  possible.  The  maximum  level  is  5 
here.  The  convergence  is  0(A *)  for  the  maximum  error  with  the  coarse  mesh 
size  going  from  No  =  80  to  160.  But  the  grid  is  approaching  a  uniform  mesh  as 
A  -»  0  and  this  slows  down  the  convergence.  As  A  -»  0  the  number  of  levels  used 
approaches  1. 

Using  method  (b),  we  let  A  -»  0  and  let  <5  =  0(A3).  Here  we  see  the  conver¬ 
gence  is  faster,  and  the  maximum  error  is  finally  0{hz).  The  ig  error  does  not 
behave  as  well. 

8.6  Estimating  the  Local  Truncation  Error  in  the  Interior.  In  [5]  we  com¬ 
pared  the  Richardson  method  and  differences  for  estimating  the  error  in  the 
interior  of  refinements  on  problem  PI.  We  found  that  there  was  very  little 
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difference  in  efficiency  between  these  methods.  The  use  of  differences  was 
slightly  more  efficient.  But  the  greater  convenience  of  Richardson  for  interior 
approximations  far  outweighs  any  small  efficiency  differences. 

B.7  Estimating  the  Local  Truncation  Error  at  Boundaries.  In  this  section  we 
will  vary  our  boundary  approximation  and  our  method  of  error  estimation  at  the 
boundary,  while  using  the  Richardson  method  in  the  interior  on  problem  P2. 

We  will  use  upwind /downwind  differencing  and  first-order  extrapolation. 
For  the  former  we  will  estimate  the  error  by  using  both  the  modified  Richardson 
method  (Section  3.3)  and  differences.  For  extrapolation  we  can  only  use 
differences.  The  results  are  shown  in  Table  3.  In  all  cases,  the  number  of  inter¬ 
vals  on  the  coarsest  mesh  is  80.  the  maximum  number  of  refinement  levels  is  5, 
and  the  refinement  ratios  N  -  M  -  4.  In  all  cases  the  fifth  refinement  level  was 
not  used.  U/D  signifies  upwind /downwind  differencing,  and  Rich,  signifies  the 
modified  Richardson  method.  The  memory  occupied  by  solution  values  is  the 
maximum  total  over  all  refinement  levels  for  one  component  of  the  solution.  All 
other  parameters  not  shown  are  the  same  as  in  the  computation  for  problem  P2 
in  Table  1. 

Clearly,  the  different  boundary  approximations  and  error  estimation 
methods  produce  approximately  the  same  results.  This  supports  our  claim  that 
our  method  of  adaptively  handling  boundaries  is  quite  general. 

6.8  How  Often  Should  the  Local  Truncation  Error  Be  Checked?.  In  Section 
2  we  used  subsequences  to  describe  the  times  at  which  we  estimate  the  local 
truncation  error  (and  possibly  alter  refinements).  In  this  section  we  shall  show 
that  for  Problem  Pi  it  is  unwise  to  monitor  the  local  truncation  error  more  often 
than  every  coarse  time  step. 
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Table  4  shows  the  results  of  these  computations  for  Problem  PI.  All  param¬ 
eters  not  mentioned  are  the  same  as  in  the  computations  for  Table  1.  The 
meaning  of  (a)  under  "tolerance  frequency”  in  Table  4  is  how  many  coarse  time 
steps  occur  between  checks  of  the  local  error.  The  column  (b)  has  two  different 
meanings,  depending  on  column  (a).  If  column  (a)  is  1  then  we  check  the  trun¬ 
cation  error  at  any  time  a  refinement  whose  level  is  less  than  or  equal  to  (b)  is 
about  to  be  advanced.  Thus,  in  these  cases  we  check  more  often  then  every 
coarse  time  step.  Table  4  shows  that  this  is  very  costly  and  produces  no  benefits 
whatever. 

If  (a)  in  Table  4  is  greater  than  one,  a  one  in  column  (b)  signifies  that  we 
check  ail  refinements  every  (a)  coarse  time  steps.  If  column  (a)  is  greater  than 
one  and  (b)  is  greater  than  one,  we  check  refinements  with  levels  greater  than 
or  equal  to  (b)  every  coarse  time  step,  and  all  others  every  (a)  coarse  time 
steps.  Of  course,  in  all  cases  in  this  table,  the  buffers  mentioned  in  Section  2.6 
have  to  be  modified,  in  a  way  analogous  to  the  argument  given  there. 

Our  results  for  these  cases  show  very  little  difference  from  checking  every 
coarse  time  step,  until  the  checking  frequency  becomes  too  seldom  (as  in  case 
(a)  =  6.  (b)  =  1).  Then  the  accuracy  starts  to  deteriorate,  because  a  pulse  may 
enter  the  boundary  before  it  is  enclosed  in  refinement(s).  (The  algorithm  could 
easily  be  modified  to  check  the  boundaries  at  every  coarse  time  step,  but  we  did 
not  do  this.) 

We  conclude  that  for  this  problem  we  may  as  well  check  the  local  error 
every  coarse  time  step,  although  this  may  depend  on  factors  such  as  the  spacing 
of  the  coarse  mesh,  the  wave  speed,  and  the  presence  of  forcing  functions 
(terms  kF  in  (2.1)).  Also,  the  results  may  be  radically  different  in  more  than 
one  space  dimension  [2], 


8.0  linear  vs.  Quadratic  Interpolation.  An  implementation  detail  we  con¬ 
sidered  is  whether  to  use  linear  or  quadratic  interpolation  when  a  level  l 
refinement  moves  into  a  region  formerly  occupied  only  by  a  level  1-1 
refinement.  In  [5]  we  found  that  there  is  practically  no  difference.  We  used  qua¬ 
dratic  interpolation  elsewhere  in  this  section. 
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Table  3.  Error  Estimation  at  Boundaries 
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Table  4  How  Often  Should  the  Local  Truncation  Error  Be  Checked? 


UNREFINED  -  REFINED 

SECOND  ORDER  UAUE  EQUATION,  T  *0 


».M  l.M  l.N  8.M  8.5* 

UNREFINED  -  REFINED 

SECOND  ORDER  UAUE  EQUATION#  T  -0 


tIM 


UNREFINED  -  REFINED  - 

SECOND  ORDER  UAVE  EQUATION,  T  *1.00 

Figure  5(c) 
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Figure  5(j) 
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UNREFINED  -  REFINED 
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